# Directed percolation in non-unitary quantum cellular automata

Ramil Nigmatullin<EMAIL_ADDRESS>Elisabeth Wagner  Gavin K.
Brennen  Center for Engineered Quantum Systems, Dept. of Physics & Astronomy,
Macquarie University, 2109 NSW, Australia

###### Abstract

Probabilistic cellular automata (CA) provides a classic framework for studying
non-equilibrium statistical physics on a lattices. A notable example is the
Domany-Kinzel CA, which has been used to investigate the process of directed
percolation and the critical dynamics of the non-equilibrium phase transition
betweeen absorbing and percolating phases. In this work, we construct a non-
unitary Quantum Cellular Automaton that generalises the Domany-Kinzel cellular
automaton and study the resulting dynamical evolution using the numerical
simulations using the tensor network iTEBD algorithm. We demonstrate the
system undergoes the absorbing/percolating phase transition and the addition
of the Hamiltonian generates coherences, which are a distinct feature of the
quantum dynamics. A proposal for the implementation of the model with Rydberg
array is put forward, which does not require local addressing of individual
sites.

##  I Introduction

In recent years there have been great advances in the development of quantum
simulation platforms. These include ultracold atoms, ions, superconducting
qubits and photonic systems. One of the most recent advances has been the
development of a quantum simulators based on arrays of ultracold Rydberg atoms
[1]. Atoms that are excited to Rydberg states interact strongly, which makes
it possible to engineer controlled multiqubit dynamics. Rydberg arrays provide
excellent platform for studying non-equilibrium dynamics of quantum many-body
systems, for example, for probing non-equilibrium quench dynamics in Ising
spin chains [2, 3], energy transport [4, 5] and the validity of thermalization
hypothesis [6, 7].

In the domain of non-equilibrium physics, one promising application of Rydberg
arrays is in exploration of non-equilibrium phase transitions (NEPTs), the
study of which is particularly challenging due to the requirement of large
system sizes and long time evolution needed to reach steady states. The
classical NEPTs are classified into universality classes characterized by
power scaling laws at critical points [8, 9]. Quantum fluctuations can
drastically alter the nature of NEPT, for example, it was shown that the
critical exponents of NEPT in a quantum analogue of the contact process model
are different to the critical exponents of a classical model [10]. The quantum
generalization of NEPT theory remains largely unexplored both theoretically
and experimentally, because of the challenges in simulating the dynamics of
quantum many-body system and engineering such systems in the lab.

Perhaps the most fundamental NEPT is the phase transition between absorbing
and percolating phase in lattice systems undergoing directed percolation (DP)
dynamics [8]. A model generalizing DP into the quantum domain, which can be
implemented using 2d Rydberg arrays, has been recently proposed [8] and
studied numerically [11, 12]. The model in [8] involves application of local
3-site unitary gates, which sequentially updates the rows, such that the state
of each row depends on the state of the previous row and corresponds to the
state of a probabilistic cellular automaton undergoing DP dynamics.

In this paper, we construct a quantum generalization of DP, which uses a
combination of non-trivial dissipative dynamics as well as coherent nearest
neighbor dynamics. The scheme is inspired by the recently developed techniques
for engineering dissipative interactions in the Rydberg arrays [13]. The
strength of the influence of quantum effects can be varied by tuning the
relative contribution of coherent and stochastic terms to the overall
dynamics. The proposed scheme has several advantages. In particular, it does
not require addressing of the individual lattice sites as all of the sites can
be updated simultaneously with spatially homogeneous dynamics. The scheme can
also be viewed as a robust way of preparing many-body quantum states, since
the reached steady states are independent of the initially prepared state. The
approach can be easily generalized to (2+1)-dimensional quantum cellular
automata for probing more complex non-equilibrium models.

The paper is organized as follows. Section II describes the construction of
the non-unitary QCA rules that generate the DP process. Section III describes
the simulation method and section IV presents the results and draw the
comparison with the classical stochastic DP. Finally, the possible physical
implementation in Rydberg arrays is discussed in section V.

##  II Model

We consider a 1-dimensional spin lattice, which undergoes open quantum
dynamical evolution modeled using a Markovian master equation. Both the
Hamiltonian and the Lindblad jump operator in the master equation are 3-qubit
operators coupling the nearest neighboring sites in the array. The master
equation is given by

$\partial_{t}\rho=\mathcal{L}[\rho]=-i[H,\rho]+\mathcal{D}[\rho],$ (1)

where the Hamiltonian is of the form

$H=\frac{1}{2}\sum_{j}\sum_{\alpha\beta}\theta_{\alpha\beta}\outerproduct{\alpha}{\alpha}_{j-1}\otimes
X_{j}\otimes\outerproduct{\beta}{\beta}_{j+1},$ (2)

where $\alpha,\beta\in\\{0,1\\}$, $\outerproduct{\alpha}{\alpha}_{j}$ is a
projector on site $j$, $X_{j}$ is the Pauli-X operator on site $j$ and
$\theta_{\alpha\beta}$ are constants whose values can be experimentally tuned.

The dissipator is given by

$\displaystyle\mathcal{D}[\rho]=$ $\displaystyle\sum_{j}\left[L^{+}_{j}\rho
L^{+\dagger}_{j}-\frac{1}{2}\left(L^{+\dagger}_{j}L^{+}_{j}\rho+\rho
L^{+\dagger}_{j}L^{+}_{j}\right)\right.$ (3)
$\displaystyle\left.+L^{-}_{j}\rho
L^{-\dagger}_{j}-\frac{1}{2}\left(L^{-\dagger}_{j}L^{-}_{j}\rho+\rho
L^{-\dagger}_{j}L^{-}_{j}\right)\right]$

with the jump operators

$L_{j}^{\pm}=\sum_{\alpha\beta}\sqrt{\gamma_{\alpha\beta}^{\pm}}\outerproduct{\alpha}{\alpha}_{j-1}\otimes\sigma^{\pm}_{j}\otimes\outerproduct{\beta}{\beta}_{j+1},$
(4)

where $\sigma^{\pm}_{j}=X_{j}\pm iY_{j}$ is a lowering/raising operators on
site $j$, $\gamma_{\alpha\beta}^{\pm}$ are tunable rate constants. If the
Hamiltonian term is set to zero, $\theta_{\alpha\beta}=0$ for all $\alpha$,
$\beta\in\\{0,1\\}$, then the master equation describes classical stochastic
dynamics on a lattice since no coherence can be generated if the system is
initialized in a product state. Such dynamics can be realized in atomic
systems such as trapped ions and Rydberg atom arrays. A possible
implementation in a Rydberg array will be described in section V.

Figure 1: A schematic representation of (1+1)-dimensional block partitioned CA
with nearest neighbor update rule, e.g. Domany-Kinzel CA. The state of a
system at time $t$ is given by contracting all the transfer matrix tensors
$\tau$ in the preceding time steps.

We will now address the question of how to appropriately choose the rates
$\gamma_{\alpha\beta}^{\pm}$, such that the model would correspond to the
process of directed percolation (DP) in the fully stochastic incoherent
regime. The idea is to choose the Lindblad operators such that the generated
dynamics can be identified with the evolution of the stochastic Domany-Kinzel
cellular automaton (DKCA) [14]. DKCA is a model exhibiting the NEPT in the DP
universality class, which is particularly simple to simulate computationally
and has been used extensively to study the absorbing/percolating critical
point. DKCA is a (1+1)-dimensional discrete stochastic CA, where the state of
each cell is updated with a probability that depends on the state of its left
and right neighbors. The probabilistic update rule
$p(s_{j}^{t}|s_{j-1}^{t},s_{j+1}^{t})$ is given by $p(1|0,0)=x$,
$p(1|0,1)=p(1,1,0)=y$, $p(1|1,1)=z$ and
$p(0|s_{j-1}^{t},s_{j+1}^{t})=1-p(1|s_{j-1}^{t},s_{j+1}^{t})$. DKCA is block
partitioned such that the updates of its even and odd cells is alternated, as
shown in figure 1. When $x=0$, the dynamics of the DKCA corresponds to DP. It
is not hard to see that setting $y=z=p$ maps the dynamics to (1+1)-dimensional
site percolation processes, where $p$ is the probability presence of a site.
The model corresponds to bond percolation, when setting $y=q$ and $z=q(2-q)$,
where $q$ is the probability of presence of a bond.

For our model to correspond to DKCA stochastic rules, we will require that the
site occupations in the stationary state of the CP map acting on an isolated
3-cell (a site and its neighbors) are the same as in the stationary states of
the DKCA rule. Let $T$ be the transfer matrix of a DKCA 3-cell rule. Let v be
the eigenvector of $T$ with the identity eigenvalue, $T\textbf{v}=\textbf{v}$.
The vector v is the stationary state of the DKCA 3-cell rule. The stationary
state of the CP map generated by the master equation satisfies
$e^{-\mathcal{L}t}[\rho^{s}]=\rho^{s}$ for all $t$; or, equivalently,
$\mathcal{L}[\rho^{s}]=0$. We choose parameters
$\\{\gamma_{\alpha\beta}^{\pm},\theta_{\alpha\beta}\\}$ such that the diagonal
elements of $\rho^{s}$ are equal to the elements v, i.e.
$\textrm{diag}(\rho^{s})=\textbf{v}$. The rates $\gamma^{\pm}_{\alpha\beta}$
are then adjusted according to

$\displaystyle\gamma^{+}_{\alpha\beta}=\begin{cases}\frac{p_{\alpha\beta}}{1-p_{\alpha\beta}}\,\gamma_{\alpha\beta}^{-}&\text{if
}\theta_{\alpha\beta}=0\\\
\frac{(1-2p_{\alpha\beta})\gamma^{-}_{\alpha\beta}+\sqrt{(\gamma^{-}_{\alpha\beta})^{2}-16p_{\alpha\beta}\theta^{2}_{\alpha\beta}(1-3p_{\alpha\beta}+2p_{\alpha\beta}^{2})}}{2(1-p_{\alpha\beta})}&\text{if
}\theta_{\alpha\beta}\neq 0\end{cases},$ (5)

whose derivation is presented in appendix A.

Note that for the model to be physical, the jump rates have to be positive and
real. This puts a restriction on the range of allowed values of
$p_{\alpha\beta}$ and $\theta_{\alpha\beta}$. If we set
$\gamma_{\alpha\beta}^{-}=1$, then we find that $\gamma_{\alpha\beta}^{+}$ is
real and positive for $p_{\alpha\beta}\geq\frac{1}{2}\
\forall\,\theta_{\alpha\beta}$, or for $p_{\alpha\beta}<\frac{1}{2}$ if
$\theta_{\alpha\beta}\leq\frac{1}{16(1-3p_{\alpha\beta}+2p_{\alpha\beta}^{2})}$.

Figure 2: Block partitioning scheme of the proposed non-unitary QCA model of
directed percolation. Partitions $\mathcal{P}_{1}$ and $\mathcal{P}_{4}$
update even cells and partitions $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$
update odd cells.

Having constructed the master equation whose dynamics implements the quantum
analogue of the DKCA 3-cell update rule, the remaining task is to define the
order in which the updates are to be carried out. We propose to split the
updates into four partitions, $\mathcal{P}_{1}$, $\mathcal{P}_{2}$,
$\mathcal{P}_{3}$ and $\mathcal{P}_{4}$ as shown in figure 2. All the cells in
each partition can be updated in parallel, since the 3-cell CP maps with each
partition act on distinct cells with no overlaps.

##  III Numerical method and calculation setup

The simulation method used in this paper is based on representing the state of
the system as matrix product state (MPS) implementing time evolution using
infinite time-evolving block decimation (iTEBD). MPS methods have been
demonstrated to be efficient tools for studying one dimensional quantum
systems [15]. If a system has translational invariance then the iTEBD method
can be used to compute the dynamics of effectively infinite system. Though
iTEBD has mostly been used to compute the unitary dynamics of closed quantum
system, it can also be applied to open dynamics with dissipation and
decoherence [16].

To apply iTEBD to Lindblad dynamics given by equation (1), we will vectorize
the density matrix using the Choi-isomorphism
$\outerproduct{\alpha}{\beta}\rightarrow\ket{\alpha}\otimes\ket{\beta}$.
Mapped in this way the density matrix
$\rho(t)=\sum_{ij}\rho_{ij}(t)\outerproduct{i}{j}$ becomes a vector in a
doubled-space $\ket{\rho(t)}=\sum_{ij}\rho_{ij}(t)\ket{i}\otimes\ket{j}$.
Under this mapping the Lindblad equation (1) becomes

$\frac{d}{dt}\ket{\rho(t)}=\mathbb{L}\ket{\rho(t)},$ (6)

where the Lindblad map $\mathbb{L}$ in doubled space is given by

$\mathbb{L}=\mathbb{H}+\mathbb{D},$ (7)

with

$\mathbb{H}=-i(H\otimes I-I\otimes H^{T}),$ (8)

and

$\mathbb{D}=\sum_{\pm}\sum_{j}\left(L_{j}^{\pm}\otimes(L_{j}^{\pm})^{*}-\frac{1}{2}\left[(L_{j}^{\pm})^{\dagger}L_{j}^{\pm}\otimes
I+I\otimes(L_{j}^{\pm})^{T}(L_{j}^{\pm})^{*}\right]\right).$ (9)

Solution to equation (1) is given by
$\ket{\rho(t)}=e^{\mathbb{L}t}\ket{\rho(0)}$, where $\ket{\rho(0)}$ is the
initial state in the vectorized representation. Defining $\mathbb{L}_{j}$ as
the Liouvillian acting on the cells in the partition $\mathcal{P}_{j}$, the
full round of the discrete non-unitary QCA model updating the state at step
$s$ to step $s+1$ is given by

$\ket{\rho_{s+1}}=e^{\mathbb{L}_{4}\tau}e^{\mathbb{L}_{3}\tau}e^{\mathbb{L}_{2}\tau}e^{\mathbb{L}_{1}\tau}\ket{\rho_{s}},$
(10)

where the time $\tau$ is the time for which the system continuously evolves
before switching to the next partition. If $\tau$ is large relative to the
characteristic time of the dissipation, i.e. $\tau\gg 1/\sqrt{\gamma}$, then
the three cells undergoing the dynamical evolution reach a steady state, which
corresponds to the state of the DKCA in the purely stochastic regime. If
$\tau$ is small, i.e. $\tau\rightarrow 0$, then the partitioning can be viewed
as Trotter decomposition of the time evolution $e^{\mathbb{L}t}$, where
$\mathbb{L}$ acts on all sites simultaneously. Thus, in the limit
$\tau\rightarrow 0$ the model can be regarded as a continuous non-unitary QCA,
as the evolution happens continuously rather than in discrete steps while the
notion of locality is still being preserved, as both the Hamiltonian and the
Lindblad jump operators are supported over nearest neighbors.

(a)  (b)

Figure 3: (a) Tensor network representation of the single round of the QCA
update with tensors $V$ and $W$ given by $V\equiv
e^{\mathbb{L}_{2}\tau}e^{\mathbb{L}_{1}\tau}$ and
$W=e^{\mathbb{L}_{4}\tau}e^{\mathbb{L}_{3}\tau}$. (b) MPS representation of
the state vector $\ket{\rho}$. The physical index dimension is $d=16$ and the
bond dimension $D$ is chosen to obtain the desired accuracy and speed of the
iTEBD algorithm.

Figure 3 shows the tensor network representation of the time evolution, where
the operator $V\equiv e^{\mathbb{L}_{2}\tau}e^{\mathbb{L}_{1}\tau}$ and
$W=e^{\mathbb{L}_{4}\tau}e^{\mathbb{L}_{3}\tau}$. The state $\rho$ is
represented as an MPS with dimension $D$. Application of the $V$ and $W$
completely positive maps increases the bond dimension of $\rho$. In TEBD
algorithms the growth of the bond dimension is controlled by only keeping a
specified number of Schmidt values. From the approximate $\ket{\rho(t)}$, one
can then compute the observables $\hat{O}$

$O(t)=\langle\hat{O}\rangle=\mathrm{tr}\left[\rho(t)\hat{O}_{D}\right]=\langle\mathbb{I}|\hat{O}_{D}|\rho(t)\rangle,$
(11)

where $\hat{O}_{D}=\hat{O}\otimes I$ and $\ket{\mathbb{I}}$ is the double
space representation of the identity operator.

##  IV Simulations results and discussion

###  IV.1 Purely stochastic limit

To verify that our non-unitary QCA model indeed represents the processes of
directed percolation on a lattice, we first compute the phase diagram of the
model in the purely stochastic regime, where $\theta_{\alpha\beta}=0$ for all
$\alpha$, $\beta$. The initial state is set to
$\rho=\outerproduct{1}{1}^{\otimes^{\infty}}$ i.e. fully active state. The
time evolution is computed using the iTEBD, for various decay rates
$\gamma^{+}$ and $\gamma^{-}$. We set $\gamma^{+}=\gamma^{-}=p$ and vary the
parameter $p$, which would correspond to the site percolation process. Both
continuous and discrete QCA model is considered by simulating the dynamics
with a 3-cell rule generated by Liouvillian evolution with short and long
durations. The lowest order Trotter error is proportional to the commutators
of the jump operators $L^{\pm}$ and hence to keep this error constant for the
continuous dynamic simulations across a range of $p$ values, we adjust the
time steps such that $\tau^{2}p(1-p)=C$ with the constant $C$ set to 0.0025.
For discrete dynamics simulations, we set $\tau=10.0$, which is sufficient
time for a 3-cell to reach a steady state before a switch to the subsequent
partitioning.

(a)

(b)

(c)

Figure 4: (a) Average number density $n$ (a)-(b) as a function of time $t$ in
the small time limit, continuous dynamics limit,
$\tau=\sqrt{0.0025/(p(1-p))}$, and no coherent dynamics, $\Omega=0$. (b) The
dynamic critical exponent $\delta$ estimated to be $0.12$ by fitting
$n(t)\propto e^{-\delta t}$ at the critical point. (c) Entropy $S(t)$ as a
function of time. The system was evolved using iTEBD with bond dimension
$D=256$ from the fully occupied initial state.

(a)

(b)

(c)

Figure 5: (a) Average number density $n$ (a)-(b) as a function of time $t$ in
the continuous dynamics limit, $\tau=10.0$, and no coherent dynamics,
$\Omega=0$. (b) The dynamic critical exponent $\delta$ estimated to be $0.19$
by fitting $n(t)\propto e^{-\delta t}$ at the critical point. (c) Entropy
$S(t)$ as a function of time. The system was evolved using iTEBD with bond
dimension $D=500$ from the fully occupied initial state.

From the MPS representing the state of the system at time $t$, we compute two
quantities - the average occupation density $n(t)$ and half-chain von Neumann
entropy $S$, the evolution of which is well understood for the DP processes
[17]. The occupation density, given by
$n(t)=\text{tr}\left(\hat{n}\rho\right)$ with $\hat{n}=\ket{1}\bra{1}$, can be
considered as an order parameter for the transition between absorbing and
percolating phases since in the absorbing phase it is zero whereas in
percolating phase it is finite. The half-chain entropy is given by
$S=(S_{1}+S_{2})/2$ with $S_{1,2}=-\text{tr}[\rho_{1,2}\log\rho_{1,2}]$, where
$\rho_{1}$ and $\rho_{2}$ are reduced density matrices of half-chain obtained
by partitioning at A-B and B-A MPS bonds, respectively. The entropy is
computed by first putting the MPS into the canonical form and then using
equations $S_{1}=-\sum_{j}^{D}{\lambda^{\prime
2}_{AB}}_{j}\log{\lambda^{\prime 2}_{AB}}_{j}$ and
$S_{2}=-\sum_{j}^{D}{\lambda^{\prime 2}_{BA}}_{j}\log{\lambda^{\prime
2}_{BA}}_{j}$, where
${\lambda^{\prime}_{AB}}_{j}={\lambda_{AB}}_{j}/\sum{\lambda_{AB}}_{j}^{2}$
and
${\lambda^{\prime}_{BA}}_{j}={\lambda_{BA}}_{j}/\sum{\lambda_{BA}}_{j}^{2}$.
The evolution of the half-chain entropy is qualitatively different in the
absorbing and percolating regimes. The entropy is zero in the absorbing steady
state since the absorbing state is a separable product state
$\rho=\ket{0}\bra{0}^{\otimes^{\infty}}$, while in the percolating phase, the
entropy of the steady state is finite reflecting the presence of classical
correlations.

Figure 4 shows the results of the simulations for the continuous QCA model.
The results indicate that there is a qualitative change in evolution at $p\sim
0.7$. In the long-time limit the occupation density tends to zero for
$p\lesssim 0.7$ and to a finite value for $p\gtrsim 0.7$ as can be seen in
figure 4 (a). The QCA evolution is evaluated for a finite number of time-steps
and the simulation does not always reach a steady state in the simulated time,
particularly, for values of $p$ close to the critical point where the
equilibration time is very long. Nevertheless, the transition between
absorbing and percolating phase is still clearly detected by evaluating the
sign of the curvature of $n(t)$ on a log-log scale i.e. for large $t$, $d\log
n/dt<0$ in the absorbing phase, $d\log n/dt>0$ in the percolating phase and
$d\log n/dt=0$ at criticality. In figure 4 (c), one can see that in the
absorbing phase the entropy $S$ initially increases as the correlations are
built up in the system, but then falls off to zero as the system approaches
the absorbing steady state. The entropy increases monotonically from zero to a
finite value in the percolating phase.

Many non-equilibrium models, including DP, can be described using the
phenomenological scaling theory. At the critical point, macroscopic
observables, such as correlation length and occupation density, are determined
by the power scaling laws with characteristic critical exponents. To determine
how close the dynamics of our QCA model is to the DP universality class, we
have computed the dynamic critical exponent $\delta$, which characterized the
critical slowing down of $n(t)$ i.e. $n(t)\propto t^{-\delta}$ at $p=p_{c}$.
The exponent $\delta$ is computed by first locating the critical point and
then doing a linear fit of the $\log n(t)$ vs $t$ curve at long time $t$ where
the slope of the curve is linear. The set of $n(t)$ curves in the vicinity of
the critical point for continuous dynamics is shown in figure 4 (b). From this
simulation data, we estimate $p_{c}=0.695$ and $\delta=0.12$. This is close
but not identical to the known values for the site DP, which are
$p^{DP}_{c}=0.705$ and $\delta^{DP}=0.16$ [8].

The results of the simulations of the dynamics in the discrete QCA limit of
large $\tau$ are shown in figure 5. The evolution of the average occupation
density in the discrete model has the same dynamical features as the
continuous model, namely, $n(t)$ decays to zero for $p<p_{c}$ and to a finite
value for $p>p_{c}$ with $p_{c}\sim 0.7$. The half-chain entropy has also a
characteristic rising and falling behavior in the absorbing phase, and
monotonically rising behavior in the percolating phase as shown in figure 5
(c). Since the Liouvillian acts for a longer time during each update step, it
takes fewer steps to reach steady state and the state update is clearly
discrete and discontinuous. We find that a higher bond dimension needs to be
used for discrete model iTEBD simulation, in order to achieve the same
simulation accuracy in the continuous case. For the discrete model, we find
that $p_{c}=0.71$ and the $\delta=0.19$, which is close to $p^{DP}_{c}=0.705$
and $\delta^{DP}=0.16$.

###  IV.2 Adding the Hamiltonian

(a)

(b)

(c)

Figure 6: Comparison of the evolution of (a) number density, (b) entropy and
(c) single site coherence, $|\rho_{01}^{s}|+|\rho_{10}^{s}|$, for different
$p$ in the presence ($\Omega=0.1$) and absence ($\Omega=0$) of the Hamiltonian
term. The results are for $\tau=\sqrt{0.0025/(p(1-p))}$; the curves $n(t)$,
$S(t)$ for $\Omega=0$ are shown in figure 4.

Having established that our non-unitary QCA model indeed corresponds to
directed percolation on a lattice, we now consider adding a local Hamiltonian
term to the Liouvillian that generates the dynamical rule. The Hamiltonian can
produce coherences and entanglement in the system, and therefore by turning it
on, one can systematically explore the influence of quantum effects on the DP
non-equilibrium phase transition. A Hamiltonian term of the form
$\outerproduct{0}{0}\otimes X\otimes\outerproduct{0}{0}$ is not compatible
with the DK model of DP since the rate $\gamma_{00}$ given by equation (5)
cannot be made positive with $p_{0,0}=0$. The local Hamiltonian term,
$I\otimes X\otimes I$, is also incompatible with DP since its basis
decomposition includes term $\outerproduct{0}{0}\otimes
X\otimes\outerproduct{0}{0}$. For this reason, we will explore the effect of
the Hamiltonian, which acts on a central qubit conditioned on the neighbors
being in state state $\ket{1}$. The Hamiltonian, $H$, is given by

$H=\Omega\outerproduct{1}{1}\otimes X\otimes\outerproduct{1}{1},$ (12)

which corresponds to equation (2) with $\theta_{11}/2=\Omega$ and all other
$\theta_{\alpha\beta}$ are set to zero.

The jump rate constants are set to $\gamma_{00}^{-}=1$, $\gamma_{00}^{+}=0$,
$\gamma_{01}^{-}=\gamma_{10}^{-}=1-y$, $\gamma_{01}^{+}=\gamma_{10}^{+}=y$,
and $\gamma_{11}^{-}=1-z$, and
$\gamma_{11}^{+}=\frac{(2z-1)\gamma_{11}^{-}+\sqrt{(\gamma_{11}^{-})^{2}-16\Omega^{2}(1-3z+2z^{2})}}{2(1-z)}$.
With these rates, the site occupations of the stationary state of the 3-cell
rule are identical to the classical DKCA of site DP, where $x=0$, $y=z=p$.

The results of the iTEBD simulation for the model in the small $\tau$ limit
with the Hamiltonian, $\Omega=0.1$, and no Hamiltonian, $\Omega=0$, for
various values of $p$ are shown in figure 6. As can be seen from figure 6 (a)
and (b), the evolution of the number densities and the half-chain entropy in
the presence of the Hamiltonian are nearly identical, with small differences
of the order of $10^{-3}$. The effect of the Hamiltonian dynamics is, however,
manifested much more evidently by the presence of coherences in the state of
the system. As a measure of coherence in the system, we will use the
${{\ell}}_{1}$-norm of coherence, given as
${C}_{1}(\rho)={\min}_{\sigma\in{\mathcal{I}}}{\parallel\rho-\sigma\parallel}_{{{\ell}}_{1}}={\sum}_{j\neq
k}|{\rho}_{j,k}|$, where ${\mathcal{I}}$ is the set of all incoherent states
and ${\parallel\cdot\parallel}_{{{\ell}}_{1}}$ is the ${{\ell}}_{1}$ matrix
norm [18]. Figure 6 (c) shows the time evolution of the single site coherence
$C_{1}$ for several values of $p$. As expected, for $\Omega=0$ the coherence
is always zero. For $\Omega=0.1$, the coherence is finite during the evolution
and its dynamics depends on whether the system is in absorbing, $p<0.7$, or
percolating phase, $p>0.7$. For $p<0.7$, the coherence rises initially and
then decays towards zero. This is because the Hamiltonian does not act on the
absorbing state $\rho=\outerproduct{0}{0}^{\otimes^{\infty}}$, and, as the
system approaches the absorbing state, the dissipation causes the coherences
to gradually decay. For $p>0.7$, the coherence rises and remains finite in a
percolating steady state. In the percolating steady state, the rate of
generation of coherence and its rate of decay are dynamically balanced.

Apart from the coherence, an important quantum feature in a many-body
dynamical system is entanglement. We have investigated the possibility of
entanglement generation by the Hamiltonian given by (12) by computing
concurrence [19, 20] in a reduced density matrix of the two adjacent sites,
which is a good measure of entanglement for mixed quantum states. We found
that the application of a single non-unitary DK QCA 3-cell rule on a product
state
$\outerproduct{1}{1}\otimes\outerproduct{1}{1}\otimes\outerproduct{1}{1}$
results in a state with non-zero concurrence. However, the concurrence
disappears after applying the gates in sequence over the four partitions
$\mathcal{P}_{1}$, $\mathcal{P}_{2}$, $\mathcal{P}_{3}$ and $\mathcal{P}_{4}$,
and there is no many-body entanglement present in the steady state in our non-
unitary QCA model.

##  V Physical implementation

The non-unitary QCA model can be effectively realized experimentally in a one
dimensional array of Rydberg atoms or ensembles. In particular, to implement
our non-unitary QCA model of DP one can adapt the set-up [13]. In [13] each
site of the one dimensional array is a 3-level system consisting of a ground
state $\ket{g}$, strongly interacting Rydberg state $\ket{r}$ and an excited
state $\ket{e}$. The sites are equidistant and the interactions are restricted
to nearest neighbors. Two laser fields couple the
$\ket{g}\leftrightarrow\ket{r}$ and the $\ket{r}\leftrightarrow\ket{e}$
transitions. The state $\ket{g}$ and $\ket{r}$ will correspond, respectively,
to the empty and active states of the DP model. The energy of the Rydberg
state depends on the state of the nearest neighbors and the energy level
$\ket{r}$ is split into three sublevels $\ket{r_{00}}$, $\ket{r_{10}}$ and
$\ket{r_{11}}$, corresponding to the state when the neighboring sites are in
states $\ket{g}\ket{g}$, $\ket{r}\ket{g}$ and $\ket{r}\ket{r}$. A
multifrequency laser field couples each of the $\ket{g}\leftrightarrow\ket{r}$
and $\ket{e}\leftrightarrow\ket{r}$ transitions with tunable coupling
strengths. The short-lived $\ket{e}$ decays to $\ket{g}$ with rate $\Gamma$
and is used to implement Lindblad jump operators, which depend on the state of
the neighbors. After adiabatically eliminating state $\ket{e}$ the effective
master equation, is defined by the Hamiltonian of the form of equation (2) and
the Lindblad operators are of the form of equation (4), except there are no
$L^{+}$ jump operator as there are no spontaneous excitation processes in the
system. The missing jump operators,
$L_{j}^{+}=\sum_{\alpha\beta}\sqrt{\gamma_{\alpha\beta}^{+}}\outerproduct{\alpha}{\alpha}_{j-1}\otimes\sigma_{j}^{+}\otimes\outerproduct{\beta}{\beta}_{j+1}$,
are needed to implement the incoherent excitation of the site $j$ conditional
on the states of the neighbors. Since $\sigma^{+}=X\sigma^{-}X$, the $L^{+}$
jump operators can be obtained by conjugating the system with the spin flip
operators and then adjusting $\gamma^{\pm}_{\alpha\beta}$ such that they
correspond to the desired jump rates. If the time between the application of
the flip pulses is small, then the action of $L^{+}$ and $L^{-}$ operators in
the master equation can be considered to be simultaneous by the Trotterization
argument.

A Rydberg atom implementation of the QCA will necessary involve a finite
number of sites. In finite size systems there is no true absorbing/percolating
phase transition as the system will always evolve towards an absorbing state
as $t\rightarrow\infty$. Nevertheless, one can still see the characteristics
of the directed percolation on the timescales which are long compared to the
duration of a single QCA step.

##  VI Conclusions

In this paper, we have constructed a one dimensional lattice model with local
unitary and non-unitary interactions, which generalizes the classical non-
equilibrium process of directed percolation into the quantum domain. The model
can be viewed as a non-unitary Quantum Cellular Automaton, with a tunable
relative strength of the coherent and incoherent dissipative dynamics. The
3-cell rule of such non-unitary QCA is the CP map generated by the open non-
unitary quantum dynamics specified by the choice of the local Lindblad jump
operators and the Hamiltonian. Varying the duration $\tau$ of the non-unitary
evolution at each time step of the QCA changes the 3-cell CP rule. We have
considered the dynamics of the model in the large $\tau$ limit, which
corresponds to discrete block partitioned non-unitary QCA, and the small
$\tau$ limit, which corresponds to continuous non-unitary QCA. The dynamics of
the model in the two limits with no Hamiltonian term was simulated using MPS
iTEBD method, and it was found that in both cases the system can undergo
absorbing/percolating non-equilibrium phase transitions. The transition point
was found to be at $p=0.71$ and $p=0.695$ for the large and small $\tau$
limits, respectively, which is within $0.7\,\%$ of the known critical value of
$p$ for the site directed percolation in (1+1)-dimensions. The dynamical
exponents $\delta$ were found to be $\sim 0.19$ and $\sim 0.12$, which appears
to be significantly different to the known critical exponent for site DP
$\delta=0.16$. However, the discrepancy could be due to the numerical error
which increases with time in the iTEBD simulations. Additional simulations
using different algorithms would be needed to obtain a more precise value of
the critical point and the exponents.

When adding a local 3-site Hamiltonian term, the dissipation rates can be
adjusted such that the evolution of the local site populations is
indistinguishable from the purely dissipative dynamics. Nevertheless, such
dynamics generates quantum coherences in the system as witnessed by the non-
zero off-diagonal elements of the density matrix. No entanglement was
generated by the chosen Hamiltonian and, thus, an interesting open question
would to determine whether certain choices of Hamiltonians and Lindblad jump
operators would generate a non-unitary QCA rules that produce entangled steady
states.

The proposed non-unitary QCA models can be engineered in ultracold atom
experiments. We have suggested a possible way of designing the rules in the
Rydberg atomic array, based on the proposal in [13]. The model in the limit of
small $\tau$ is equivalent to the continuous non-unitary QCA and, thus, can be
implemented by acting on all sites with identical laser fields. This global
control features makes the scheme attractive for experimental implementation.

Since Rydberg arrays can be configured into two dimensional lattice
configurations [1], a natural extension of the present work is to consider
(2+1)-dimensional non-unitary QCA. While the generalization of the
construction of the non-unitary QCA rules from one to two dimensions is
straightforward, the dynamics of the two dimensional CA is usually
significantly more complex [21]. The quantum mechanical generalization of two
dimensional CA models is challenging to simulate on classical computers and
for this reason their exploration are likely to be driven by experiments. Such
models would constitute an excellent test bed for existing and novel tensor
network algorithms.

##  VII Acknowledgements

This work was supported by the Australian Research Council Centre of
Excellence for Engineered Quantum Systems (Grant No. CE 170100009). This
research was undertaken with the assistance of resources from the National
Computational Infrastructure (NCI), which is supported by the Australian
Government.

## References

  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nature Physics 16, 132 (2020).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Guardado-Sanchez  _et al._ [2018] E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Probing the quench dynamics of antiferromagnetic correlations in a 2d quantum ising spin system, Phys. Rev. X 8, 021069 (2018).
  * Barredo  _et al._ [2015] D. Barredo, H. Labuhn, S. Ravets, T. Lahaye, A. Browaeys, and C. S. Adams, Coherent excitation transfer in a spin chain of three rydberg atoms, Phys. Rev. Lett. 114, 113002 (2015).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with rydberg atoms, Science 365, 775 (2019), science.sciencemag.org/content/365/6455/775.full.pdf .
  * Kim  _et al._ [2018] H. Kim, Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Detailed balance of thermalization dynamics in rydberg-atom quantum simulators, Phys. Rev. Lett. 120, 180502 (2018).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018).
  * Hinrichsen [2000] H. Hinrichsen, Non-equilibrium critical phenomena and phase transitions into absorbing states, Advances in Physics 49, 815 (2000), https://doi.org/10.1080/00018730050198152 .
  * Ódor [2004] G. Ódor, Universality classes in nonequilibrium lattice systems, Rev. Mod. Phys. 76, 663 (2004).
  * Gillman  _et al._ [2019] E. Gillman, F. Carollo, and I. Lesanovsky, Numerical simulation of critical dissipative non-equilibrium quantum systems with an absorbing state, New Journal of Physics 21, 093064 (2019).
  * Gillman  _et al._ [2020] E. Gillman, F. Carollo, and I. Lesanovsky, Nonequilibrium phase transitions in ($1+1$)-dimensional quantum cellular automata with controllable quantum correlations, Phys. Rev. Lett. 125, 100403 (2020).
  * [12] E. Gillman, F. Carollo, and I. Lesanovsky, Numerical simulation of critical quantum dynamics without finite size effects, arXiv:2010.10954 .
  * Wintermantel  _et al._ [2020] T. M. Wintermantel, Y. Wang, G. Lochead, S. Shevate, G. K. Brennen, and S. Whitlock, Unitary and nonunitary quantum cellular automata with rydberg arrays, Phys. Rev. Lett. 124, 070503 (2020).
  * Domany and Kinzel [1984] E. Domany and W. Kinzel, Equivalence of cellular automata to ising models and directed percolation, Phys. Rev. Lett. 53, 311 (1984).
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Orús and Vidal [2008] R. Orús and G. Vidal, Infinite time-evolving block decimation algorithm beyond unitary evolution, Phys. Rev. B 78, 155117 (2008).
  * Harada and Kawashima [2019] K. Harada and N. Kawashima, Entropy governed by the absorbing state of directed percolation, Phys. Rev. Lett. 123, 090601 (2019).
  * Carmeli  _et al._ [2018] C. Carmeli, T. Heinosaari, S. Maniscalco, J. Schultz, and A. Toigo, Determining quantum coherence with minimal resources, New Journal of Physics 20, 063038 (2018).
  * Wootters [1998] W. K. Wootters, Entanglement of formation of an arbitrary state of two qubits, Phys. Rev. Lett. 80, 2245 (1998).
  * Wong and Christensen [2001] A. Wong and N. Christensen, Potential multiparticle entanglement measure, Phys. Rev. A 63, 044301 (2001).
  * Rujàn [1987] P. Rujàn, Cellular automata and statistical mechanical models, Journal of Statistical Physics 49, 139 (1987).

##  Appendix A Derivation of rate constants $\gamma_{\alpha\beta}^{\pm}$

We present a derivation of the continuous-time version of the DKCA dynamics
using the Lindblad formalism. The principal idea consists of calculating the
stationary states after each time step, and setting them in such a way that
they obey the dynamics of the desired discrete system. More specifically, the
neighborhood-dependent amplitudes of the Hamiltonian ($\theta_{\alpha\beta}$)
and jump operators ($\gamma_{\alpha\beta}^{\pm}$) are derived in accordance
with the update probabilities ($p_{\alpha\beta}$) from the DKCA.

The stationary states are determined by

$\displaystyle\mathbb{L}\ket{\rho^{s}}=0,$ (13)

with $\mathbb{L}$ being the vectorized version of the Liouvillian in equations
(7) to (9) in the main text; writing

$\displaystyle\mathbb{L}$
$\displaystyle=\sum_{\alpha\beta\alpha^{\prime}\beta^{\prime}}\outerproduct{\alpha}{\alpha}_{j-1}\otimes\outerproduct{\alpha^{\prime}}{\alpha^{\prime}}_{j^{\prime}-1}$
$\displaystyle\qquad\qquad\otimes\Big{[}-i\frac{\theta_{\alpha\beta}}{2}\left(\outerproduct{0}{1}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}+\outerproduct{0}{1}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}+\outerproduct{1}{0}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}+\outerproduct{1}{0}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}\right)$
$\displaystyle\qquad\qquad\qquad+i\frac{\theta_{\alpha^{\prime}\beta^{\prime}}}{2}\left(\outerproduct{0}{0}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}+\outerproduct{1}{1}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}+\outerproduct{0}{0}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}+\outerproduct{1}{1}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}\right)$
$\displaystyle\qquad\qquad\qquad+\sqrt{\gamma_{\alpha\beta}^{+}\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}\outerproduct{0}{1}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}+\sqrt{\gamma_{\alpha\beta}^{-}\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}\outerproduct{1}{0}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}$
$\displaystyle\qquad\qquad\qquad-\frac{\gamma_{\alpha\beta}^{-}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}{2}\outerproduct{0}{0}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}-\frac{\gamma_{\alpha\beta}^{-}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}{2}\outerproduct{0}{0}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}$
$\displaystyle\qquad\qquad\qquad-\frac{\gamma_{\alpha\beta}^{+}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}{2}\outerproduct{1}{1}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}-\frac{\gamma_{\alpha\beta}^{+}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}{2}\outerproduct{1}{1}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}\Big{]}$
$\displaystyle\qquad\qquad\otimes\outerproduct{\beta}{\beta}_{j-1}\otimes\outerproduct{\beta^{\prime}}{\beta^{\prime}}_{j^{\prime}-1}$
(14)

by considering the Hamiltonian (2) and jump operators (4)
$\forall\alpha,\alpha^{\prime},\beta,\beta^{\prime}\in\\{0,1\\}$, where the
indices $j^{\prime}$ denote the position in the dual space arising from the
vectorization approach. The steady states of the system are then defined by
the system of equations

$\displaystyle\text{I}.$ $\displaystyle\qquad
0=\frac{-i}{2}(\theta_{\alpha\beta}\
\rho_{\alpha,\alpha^{\prime},1,0,\beta,\beta^{\prime}}^{s}-\theta_{\alpha^{\prime}\beta^{\prime}}\
\rho_{\alpha,\alpha^{\prime},0,1,\beta,\beta^{\prime}}^{s})-\frac{\gamma_{\alpha\beta}^{-}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}{2}\
\rho_{\alpha,\alpha^{\prime},0,0,\beta,\beta^{\prime}}^{s}+\sqrt{\gamma_{\alpha\beta}^{+}\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}\
\rho_{\alpha,\alpha^{\prime},1,1,\beta,\beta^{\prime}}^{s}$
$\displaystyle\text{II}.$ $\displaystyle\qquad
0=\frac{-i}{2}(\theta_{\alpha\beta}\
\rho_{\alpha,\alpha^{\prime},1,1,\beta,\beta^{\prime}}^{s}-\theta_{\alpha^{\prime}\beta^{\prime}}\
\rho_{\alpha,\alpha^{\prime},0,0,\beta,\beta^{\prime}}^{s})-\frac{\gamma_{\alpha\beta}^{-}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}{2}\
\rho_{\alpha,\alpha^{\prime},0,1,\beta,\beta^{\prime}}\textbf{}$
$\displaystyle\text{III}.$ $\displaystyle\qquad
0=\frac{-i}{2}(\theta_{\alpha\beta}\
\rho_{\alpha,\alpha^{\prime},0,0,\beta,\beta^{\prime}}-\theta_{\alpha^{\prime}\beta^{\prime}}\
\rho_{\alpha,\alpha^{\prime},1,1,\beta,\beta^{\prime}})-\frac{\gamma_{\alpha\beta}^{+}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}{2}\
\rho_{\alpha,\alpha^{\prime},1,0,\beta,\beta^{\prime}}^{s}$
$\displaystyle\text{III}.$ $\displaystyle\qquad
0=\frac{-i}{2}(\theta_{\alpha\beta}\
\rho_{\alpha,\alpha^{\prime},0,1,\beta,\beta^{\prime}}^{s}-\theta_{\alpha^{\prime}\beta^{\prime}}\
\rho_{\alpha,\alpha^{\prime},1,0,\beta,\beta^{\prime}}^{s})-\frac{\gamma_{\alpha\beta}^{+}+\gamma_{\alpha^{\prime}\beta^{\prime}}^{+}}{2}\
\rho_{\alpha,\alpha^{\prime},1,1,\beta,\beta^{\prime}}^{s}+\sqrt{\gamma_{\alpha\beta}^{-}\gamma_{\alpha^{\prime}\beta^{\prime}}^{-}}\
\rho_{\alpha,\alpha^{\prime},0,0,\beta,\beta^{\prime}}^{s},$ (15)

where the first and last two subscripts of the steady state elements
$\rho_{\alpha,\alpha^{\prime},c,c^{\prime},\beta,\beta^{\prime}}^{s}$ denote
the basis of the left or right sites, $(j-1,j^{\prime}-1)$ or
$(j+1,j^{\prime}+1)$ respectively, while the center site $(j,j^{\prime})$ is
associated with $c,c^{\prime}\in\\{0,1\\}$.

W.l.o.g. one can consider that the nearest-neighbor states and their dual
space correspondence are identical, setting $\alpha=\alpha^{\prime}$ and
$\beta=\beta^{\prime}$. The Liouvillian is then simplified to

$\displaystyle\mathbb{L}$
$\displaystyle=\sum_{\alpha\beta}\outerproduct{\alpha}{\alpha}_{j-1}\otimes\outerproduct{\alpha}{\alpha}_{j^{\prime}-1}$
$\displaystyle\qquad\qquad\otimes\Big{[}\frac{-i\theta_{\alpha\beta}}{2}\Big{(}\outerproduct{0}{1}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}+\outerproduct{0}{1}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}+\outerproduct{1}{0}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}}+\outerproduct{1}{0}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\outerproduct{0}{0}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}-\outerproduct{1}{1}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}-\outerproduct{0}{0}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}-\outerproduct{1}{1}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}\Big{)}$
$\displaystyle\qquad\qquad\qquad+\gamma_{\alpha\beta}^{+}\,(\outerproduct{0}{1}_{j}\otimes\outerproduct{0}{1}_{j^{\prime}}-\outerproduct{1}{1}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}})$
$\displaystyle\qquad\qquad\qquad+\gamma_{\alpha\beta}^{-}\,(\outerproduct{1}{0}_{j}\otimes\outerproduct{1}{0}_{j^{\prime}}-\outerproduct{0}{0}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}})$
$\displaystyle\qquad\qquad\qquad-\frac{\gamma_{\alpha\beta}^{-}+\gamma_{\alpha\beta}^{+}}{2}\,(\outerproduct{0}{0}_{j}\otimes\outerproduct{1}{1}_{j^{\prime}}+\outerproduct{1}{1}_{j}\otimes\outerproduct{0}{0}_{j^{\prime}})\Big{]}$
$\displaystyle\qquad\qquad\otimes\outerproduct{\beta}{\beta}_{j-1}\otimes\outerproduct{\beta}{\beta}_{j^{\prime}-1}.$
(16)

For clarity the subscripts $\alpha,\beta$ are omitted further by replacing
$\gamma_{\alpha\beta}^{\pm}\rightarrow\gamma^{\pm}$,
$p_{\alpha\beta}\rightarrow p$, $\theta_{\alpha\beta}/2\rightarrow\Omega$, and
$\rho_{\alpha,\alpha^{\prime},c,c^{\prime},\beta,\beta^{\prime}}^{s}=\rho_{\alpha,\alpha,c,c^{\prime},\beta,\beta}^{s}\rightarrow\rho_{c,c^{\prime}}^{s}$;
yielding a simplified system of equations that define the steady state
elements $\forall\alpha,\beta$:

$\displaystyle\text{I}.$ $\displaystyle\qquad
0=i\Omega(\rho_{01}^{s}-\rho_{10}^{s})-\gamma^{-}\rho_{00}^{s}+\gamma^{+}\rho_{11}^{s}$
$\displaystyle\text{II}.$ $\displaystyle\qquad
0=i\Omega(\rho_{00}^{s}-\rho_{11}^{s})-\frac{\gamma^{+}+\gamma^{-}}{2}\rho_{01}^{s}$
$\displaystyle\text{III}.$ $\displaystyle\qquad
0=i\Omega(\rho_{11}^{s}-\rho_{00}^{s})-\frac{\gamma^{+}+\gamma^{-}}{2}\rho_{10}^{s}$
$\displaystyle\text{IV}.$ $\displaystyle\qquad
0=i\Omega(\rho_{10}^{s}-\rho_{01}^{s})+\gamma^{-}\rho_{00}^{s}-\gamma^{+}\rho_{11}^{s},$
(17)

where it is to note that equations II. and III. demand
$\rho_{01}^{s}=(\rho_{10}^{s})^{*}$, and I. and IV. lead to
$\rho_{00}^{s}+\rho_{11}^{s}=1$ representing the trace-preserving property of
the reduced density matrix at site $j$.

Desired dynamics are thus encoded into the parameters that describe the
system, $\gamma^{\pm}$ and $\Omega$. As $\Omega$ describes the strength of the
Hamiltonian which tunes quantum coherences, this shall not be fixed, but
remain a variable to effectively increase or decrease the quantum coherences,
turning them on $(\Omega\neq 0)$ or off $(\Omega=0)$.

Here, the stationary states of the system are set after each time step to obey
the dynamics of the DKCA. In this model an arbitrary input state is taken to
$\begin{pmatrix}p&0\\\ 0&1-p\end{pmatrix}$, i.e. the system is in the
classical ‘1’ state with probability $p$, or in ‘0’ with probability $1-p$
after one time step. The corresponding quantum dynamics are described in the
Lindblad formalism by setting $\rho_{11}^{s}=1-p$. and then solving the system
of equations (17) for one of the variables that describe the system, say
$\gamma^{+}$. By introducing the Hamiltonian part of the Liouvillian quantum
coherences show up as the off-diagonal matrix elements of the density matrix:
$\begin{pmatrix}p&\rho_{01}^{s}\\\ -\rho_{01}^{s}&1-p\end{pmatrix}$ with

$\displaystyle\rho_{01}^{s}$ $\displaystyle=\begin{cases}0&\text{, if
}\Omega=0\\\
\frac{-i}{4\Omega}\left(\gamma^{-}-\sqrt{(\gamma^{-})^{2}-16\,\Omega^{2}\,(1-3p+2p^{2})}\right)&\text{,
if }\Omega\neq 0\end{cases},$ (18)

where

$\displaystyle\gamma^{+}$
$\displaystyle=\begin{cases}\frac{p}{1-p}\gamma^{-}&\text{, if }\Omega=0\\\
\frac{(2p-1)\gamma^{-}+\sqrt{(\gamma^{-})^{2}-16\,\Omega^{2}(1-3p+2p^{2})}}{2(1-p)}&\text{,
if }\Omega\neq 0\end{cases},$ (19)

which is e.g. with $\gamma^{-}=1$ only real and positive in the cases

$\displaystyle\text{I}.$ $\displaystyle\qquad\Omega=0$
$\displaystyle\text{II}.$ $\displaystyle\qquad p\geq\frac{1}{2}\
\forall\,\Omega$ $\displaystyle\text{III}.$ $\displaystyle\qquad
p<\frac{1}{2}\text{ if }\Omega\leq\frac{1}{16(1-3p+2p^{2})}.$ (20)

Note that for our simulations of the DP model the Hamiltonian is chosen to be
only turned on if both neighboring states are in the ‘1’-state, i.e. we set
$\Omega=\theta_{11}$/2 where $\alpha=\beta=1$.