Quantum and classical temporal correlations in $(1 + 1)D$ Quantum Cellular Automata

We employ $(1 + 1)$-dimensional quantum cellular automata to study the evolution of entanglement and coherence near criticality in quantum systems that display non-equilibrium steady-state phase transitions. This construction permits direct access to the entire space-time structure of the underlying non-equilibrium dynamics. It contains the full ensemble of classical trajectories and also allows for the analysis of unconventional correlations, such as entanglement in the time direction between the"present"and the"past". Close to criticality, the dynamics of these correlations - which we quantify through the second-order Renyi entropy - displays power-law behavior on its approach to stationarity. Our analysis is based on quantum generalizations of classical non-equilibrium systems: the Domany-Kinzel cellular automaton and the Bagnoli-Boccara-Rechtman model, for which we provide estimates for the critical exponents related to the classical and quantum components of the entropy. Our study shows that $(1 + 1)$-dimensional quantum cellular automata permit an intriguing perspective on the nature of classical and quantum correlations in out-of-equilibrium systems.

Introduction. Cellular automata (CA) are paradigmatic models for the study of non-equilibrium processes and phase transitions that fall outside the realm of equilibrium statistical mechanics [1][2][3][4]. They also serve as models for complex dynamical processes or even for computation [5,6]. Since their inception, considerable activity has been dedicated to generalizing the CA concept into the quantum domain [7][8][9][10] (see e.g. Ref. [11] for a review). These efforts have resulted in a host of different settings which combine CA dynamics with quantum coherent evolution. Here, we follow a route introduced in Refs. [12,13], which establishes a direct connection between classical probabilistic CA and the discrete-time dynamics of a quantum system on a (1+1)-dimensional lattice, i.e. (1+1)D quantum cellular automata (QCA). The underlying idea is that, for certain models, the dynamics of classical (diagonal) observables can be described through the Master equation of an associated classical non-equilibrium process. This establishes a well-defined classical limit, which serves as a convenient starting point for investigations into the impact of quantum correlations on non-equilibrium processes in many-body systems [12,[14][15][16][17][18][19]. A further appealing aspect is that these systems are realizable on current quantum simulation platforms, such as two-dimensional Rydberg lattice gases [20][21][22][23].
In this work we exploit that the (1 + 1)D QCA construction encodes an entire space-time dynamics within a single pure quantum state, i.e. the whole ensemble of possible histories (space-time trajectories) of a nonequilibrium process becomes accessible by measuring observables of the underlying two-dimensional lattice. This encoding allows to access and analyze new quantities, such as quantum entanglement in the time-domainbetween present and past. We demonstrate this for (1 + 1)D QCA are based on the successive application of unitary quantum gates. By acting on pairs of rows sequentially, an entangled 2D quantum state |ψ is created that encodes the entire space-time structure of the ensuing non-equilibrium dynamics. (b) An elementary cellular automaton (ECA) generates a 2D product state. On the QCA an ECA can be realized by gates that merely flip the unoccupied target sites into occupied ones. The image shows snapshots taken at three times (t1 < t2 < t3) for an evolution under rule 150, see Table I. (c) An evolution similar to a probabilistic classical CA is implemented when the gate rotates target sites into superposition states. This generates an entangled 2D quantum state. The image shows the density of occupied sites of the BBR model. (d) Unitary evolution of a 1D system in the (1 + 1)D QCA framework, for comparison. Here, the space-time structure for this process is inaccessible.
QCA related to quantum generalizations of two paradigmatic non-equilibrium models: the Domany-Kinzel cellular automaton (DKCA) [24] and the Bagnoli-Boccara-Rechtman (BBR) model [25]. Both systems possess a critical point associated with an absorbing phase transition (APT). In its vicinity the second-order Renyi entropy -which quantifies the present-past entanglement -displays power-law scaling on its approach to stationarity. While such scaling was observed previously in classical systems [26], we show that in the associated  Fig. 1(b) and (c) for the evolution of an initial seed generated by the ECA and BBR rules.
quantum models the contribution to the entropy stemming from quantum coherences also displays a powerlaw decay. The value of the associated critical exponent is such that at long times the entropy is dominated by the classical contribution. This yields not only interesting insights into the dynamical behavior of quantum and classical correlations, but in general highlights the potential offered by the QCA platform which grants access to the full space-time structure of quantum non-equilibrium processes.
(1 + 1)D QCA and link to other CA models. QCA are lattice systems with sites that may be either occupied or empty, see Fig. 1(a). The dynamics is constructed such that the occupation of sites in a rowor more generally in a d-dimensional surface -is determined by those of the row above it [12,13]. This leads to an effective time-dimension, and the corresponding models are termed (d + 1)-dimensional with d spatial dimensions perpendicular to the single "time" dimension [4]. Similar models have been studied extensively in the context of classical systems [27]. In (1 + 1)D QCA, at any time t the system is described by a pure quantum state |ψ t , which is an element of the Hilbert space H = τ,k H τ,k , where H τ,k are local Hilbert spaces of a 2D lattice indexed by (τ, k). Here we consider local Hilbert spaces that are two-dimensional with basis {|• , |• }, where n |• = 0, n |• = |• and n is the local particle number operator. We will refer to the states |• and |• as occupied and empty respectively.
The (1 + 1)D QCA evolves under the action of unitary gates G τ,k . These apply local updates to a "target" site at (τ, k), depending on the state of a set of "control" sites in row τ − 1 that form the "neighbourhood", N k , of site k, see Fig. 1(a). Since we are considering binary variables, N k is taken to be an integer labelling the possible neighbourhoods, whose binary representation (see Table  I) specifies the occupation of the sites [6]. The state |ψ t evolves in discrete time-steps as |ψ t = G t |ψ t−1 . Here the "global update" G t consists of an ordered product of G t,k , one per each site of row t. We will consider the situation in which the 2D state is initialised at t = 1 into a product state of all unoccupied sites, |• , except for the first row, τ = 1. This will have a single occupied site, |• , at the center, which we refer to as the seed initial condition. During the subsequent evolution, |ψ t then encodes the entire space-time structure of cellular automaton dynamics from this initial condition, i.e. it allows for access to the full history of trajectories, permitting the analysis of typically inaccessible (quantum) correlations between different space-time regions.
In the following, we consider gates that act on a single target site and a three-site "control" neighbourhood with the form, The symbol ⊗ separates the operators which act on the controls (placed to the left) and on the target (placed on the right): the P N k = |N k N k | are projectors onto the different "classical" configurations of the control sites [see Fig. 1(a) and Table I] and σ x = |• •| + |• •| acts on the target site. The angles α(N k ) determine by how much the target site is rotated from its initial state |• into a superposition. The probability for the target site to be occupied, given a particular configuration of the control sites, is P For the gate (1), a local update (or "rule") is specified by choosing the eight (real) values of α(N k ) or equivalently p(N k ). As we discuss in the following this includes several informative cases, listed in Table I, which directly connect to CCA with simultaneous updates. This is because for these cases the G τ,k for different k commute, such that all choices of G t are equivalent.
Choosing rotation angles such that p(N k ) = 0, 1, the gate of Eq. (1) reproduces deterministic CCA on the 2D state |ψ t , which remains in a product (unentangled) form at all times, see Fig. 1(b). For example, ECA can be realised in this setting including irreversible cases [6], which corresponds to the classical result that irreversible CCA can be embedded in higher-dimensional reversible CCA [29]. Further to that, Eq. (1), contains the case of classical PCA [see Fig. 1(c)]. Here the diagonal components of the density matrix Ξ t = |ψ t ψ t | are equal to the probability of producing the corresponding space-time configuration under a PCA dynamics [12]. However, the unitary dynamics of the (1 + 1)D QCA can generate offdiagonal terms in Ξ t , i.e., coherence. This means that the gate in Eq. (1) generalizes any desired PCA into a genuine quantum system: the (1 + 1)D QCA encodes the original classical dynamics, including all associated physics such as APTs, while also displaying uniquely quantum features. For example, the BBR model is a PCA with a three-site control neighbourhood that displays APTs [25,30]. It displays two absorbing states (the fully occupied and fully empty product states) and is totalistic, meaning that the local updates depend only on the total number of occupied sites in the neighbourhood, see Table  I. The corresponding totalistic update rule with a twosite neighbourhood, the DKCA, which displays a single absorbing state of all empty sites, can also be considered in this three-site neighbourhood setting, e.g., by making updates not depending on the control site in the middle. As for deterministic CCA, |ψ t encodes also in this case the entire space-time history, although now being an entangled state, see Fig. 1

(c).
Scaling of temporal correlations. In the following, we will focus on (1 + 1)D QCA with gates as in Eq. (1), that realize quantum generalizations of the DKCA and the BBR model. We will refer to these as QDKCA and QBBR model respectively. At any time, the QCA can be partitioned in such a way that the t-th row is singled out. Since rows with τ > t are in a product state, the entanglement between the two subsystems generated by the partition quantifies the amount of quantum correlations in the pure state |ψ t , between the "present" (sites in row t) and the "past" (rows τ < t). We measure this "present-past" entanglement through the second-order Renyi entropy, i.e. the logarithm of the purity of the reduced state ρ t = Tr τ =t Ξ t of the t-th row: Given that we are considering extensions of PCA, it is natural to expand ρ t into a diagonal part, ρ cl t , and an off-diagonal part, X t , such that ρ t = ρ cl t + X t . For classical processes all off-diagonal terms are zero and ρ t = ρ cl t [26,28]. Under this decomposition the purity of the reduced state becomes a sum of two terms: The first one is equivalent to the classical component of the purity. This can be viewed as due to the probabilistic nature of the process through which sites can be rotated into the occupied state |• . The second term is the 2 -norm of the density matrix coherence. This contribution is positive, and zero only if ρ t = ρ cl t . As such, it can only increase the purity of ρ t and is a manifestation of the quantum correlations present in the QCA. While sufficient for our purposes, we note that C 2 is not a strict measure of coherence as it violates certain monotonicity conditions [31].
For some classical PCA, it was recently found that near the critical point of an ATP the second-order Renyi entropy scales as S cl 2 ∼ t −p , with a universal exponent p = 0.632613(6) [26]. The same behavior is expected, by construction, when considering ρ cl t of the QDKCA and the QBBR model, and hence should display a power-law decay at criticality.
In what follows, we focus on the quantum Renyi entropy in Eq. (2) and show that, in the vicinity of the APT of the QDKCA and the QBBR model, it obeys a scaling form Furthermore, we observe that the coherence also follows a similar power-law decay close to criticality, which defines an additional critical exponent q coh . For the seed initial conditions considered, we establish that q coh q ent . As might then be expected, we find that the quantum entropy S 2 (t) tends to the one of the corresponding classical PCA, S cl 2 (t), for sufficiently long times. Numerical Results. In the following we present results concerning the scaling of S 2 and C 2 close to criticality for the QDKCA and the QBBR. We compute these quantities by simulating the reduced evolution of ρ t using tensor networks (TNs) and matrix product states (MPSs) [32,33]. The method employed -detailed in Ref.
[28] -relies on representing ρ t as an MPS and G as a matrix product operator (MPO), and applying standard methods for simulating MPS evolution [34,35]. Details concerning the lattice size L, MPS bond-dimension χ and other parameters related to the simulation are contained in the caption of Fig. 2. All simulations start with an initial seed placed in the center of the first time-slice [see Fig. 1 In Fig. 2(a) we show the entropy S 2 for the QDKCA with p 2 = 0.874 and the six values of p 1 indicated in the figure (solid lines). For the two lowest values of p 1 , the second-order Renyi entropy rapidly vanishes, due to the fact that the systems approaches an absorbing (product) state. In contrast, the curves with the three highest values of p 1 tend to stationary values. This demonstrates that S 2 can play the role of an order parameter, by distinguishing between the two different phases. When choosing p 2 = p 1 (2 − p 1 ) the DKCA is equivalent to socalled bond-directed-percolation and extensive studies of this classical process have determined the location of the critical point as p 2 = 0.874, p 1 = 0.645 [4]. The corresponding critical curve is shown as solid green lines in  Note that, in estimating critical exponents for APTs, determining the location of the critical point is typically a key source of error [4]. However, since the QDKCA and DKCA share the same critical point by construction, in our case this error is negligible. The relevant error sources for the QDKCA model are thus associated with finite-L, finite-χ and finite-time effects. In contrast, due to the relatively few studies on the BBR model, for the QBBR model the uncertainty on the location of the critical point provides a significant contribution to the error. Overall, in addition to larger finite-χ errors, the error associated to the estimates of the critical exponents for the QBBR model are considerably larger than those of the QDKCA. For details on the estimation of errors, and a further discussion of related issues, see Ref.
[28]. Fig. 2(a) also displays S cl 2 (dashed lines). For each p 1 , as t increases, the curves for S 2 and S cl 2 become indistinguishable on the scale shown. This means that the critical exponent found also holds for the classical model. We remark that this is different from the value obtained for S cl 2 in Ref. [26], which may be due to the fact that we are not using a homogeneous initial condition but an initial seed instead. The agreement between S 2 (t) and S cl 2 (t) suggests that C 2 becomes irrelevant compared with Tr (ρ cl t ) 2 over time [28]. This is confirmed in Fig. 2(b), where C 2 is shown for the same values of p 1 . As with S 2 , at criticality we observe a power-law behaviour in C 2 . However, we find that the timescale over which C 2 approaches a power-law is considerably larger than that of S 2 . As such, to estimate the critical exponent, we construct the time-dependent effective exponent, q coh (t) = − ln [C 2 (t)/C 2 (t/2)], as shown in the inset of Fig. 2(b). As can be seen, C 2 does indeed approach a power-law (indicated by q coh (t) approaching a constant value), and we estimate the exponent to be q coh = 2.96 ± 0.04 by averaging over the effective exponent between t = [450, 550]. Conclusions and Outlook. QCA constitute a platform that allows to realize a number of canonical CA scenarios. They can be experimentally realized on quantum simulators and encode the entire space-time information of a non-equilibrium process in a single quantum state. This permits experimental access to unusual properties, such as entanglement in the time domain. Already simple QCA, which are quantum generalizations of the classical DKCA and the BBR model, reveal intriguing features, such as power-law scaling of entanglement and coherence with time at criticality. In the future it would be interesting to focus on more intricate situations, e.g. QCA where the elementary gates do not commute, so that the order in which local updates are applied defines inequivalent global updates G t [13]. In such a setting, the updates in (1 + 1)D QCA can be considered as asynchronous updates, the impacts of which have been extensively studied in the classical case [36][37][38] but are still largely unexplored in the quantum domain.

REDUCED EVOLUTION OF THE 1D QUANTUM STATE OF ROW t
In this section we show how, starting from the dynamics of a (1 + 1)D QCA, it is possible to obtain the timeevolution of a 1D QCA. This is achieved by focusing on the reduced state of the last updated row at every discrete time.

Discrete update equation of 1D QCA
For the 2D state |ψ t we can construct the density matrix Ξ t = |ψ t ψ t |, so that the reduced state of row t is readily obtained as ρ t = Tr τ =t [Ξ t ], i.e. by tracing out all the degrees of freedom of the (1 + 1)D QCA but those in row t. Such a reduced state can then be related to the reduced state of the row t − 1 at time t − 1, ρ t−1 = Tr τ =t−1 [Ξ t−1 ], to define a discrete-time evolution equation for the 1D state ρ t .
To this end, we note that, from row t + 1 onward, the state |ψ t features all sites in the empty state and is thus in a product form. Furthermore, the global update gate G t acts non-trivially only on rows t and t − 1. As such, using the fact that |ψ t = G t |ψ t−1 , the reduced state ρ t can be related to ρ t−1 as, (S1) Here, |Ω t is a 1D product state of all empty sites describing the row t just before the update. This clearly shows that the evolution of the (1 + 1)D QCA induces a discrete time-evolution of a 1D quantum state ρ t . While the evolution of the 2D state |ψ t is unitary, so that the state remains always pure, the evolution of such a 1D QCA is in general non-unitary and the state ρ t mixed. In this sense, the unitary (1 + 1)D QCA can induce a non-unitary 1D QCA.

Relation to system-environment evolution
The evolution of the 1D state ρ t can be related to the usual evolution of a coupled system-environment under Markovian assumptions. For a system state ρ S and environment-state φ E , the system state at time t can be written as, where U t is the joint system-environment unitary.
Such an evolution can be made equivalent to the reduced evolution (S2) by choosing U t = SG t where S |nm = |mn is a SWAP gate acting on the system and environment. In that case, Taking ρ S = ρ t−1 and φ E = |Ω t Ω t | then reproduces the 1D state evolution of Eq. (S2).

Relation to unitary 1D QCA
As a simple example of an evolution where the purity of ρ t is preserved, one can consider a local update gate of the form G = U ⊗ 1, with U = exp [−iδt(h 1,2 + h 2,3 )], δt = 0.01 and h i,j = σ x i n j + n i σ x j . This is the choice that we made to produce the plot of Fig. 1(d). In this case we can write the global update gate as G t = U t−1 ⊗ 1 t , where U t is a unitary, ordered product of the U s (one per target site) which act on row t − 1.
It is important to notice that this choice of G t does not entangle the rows t and t − 1 (indeed in this example, the row t is not modified at all). As such, if SWAP gates are applied following the global update, then the evolution of ρ t is unitary as shown by the following iterative equation:

RELATIONS TO CLASSICAL MODEL
Pure state representation of classical model The classical probability distribution of a PCA can be encoded in a so-called probability vector [26] via the relation where one has m P t (m) = 1. Here, indeed, m labels the set of orthonormal states of definite occupation (i.e. the classical configurations), while P t (m) is the probability of the state m occurring in the PCA at time t. The norm of this state provides the "purity" of such a classical probability distribution, that we call γ cl , The value of γ cl is 1, denoting purity of the state, if and only if the probability is such that P t (m) = 1, only for a single state |m and zero for all the others. Equivalently, the distribution P t (m) can be encoded in a diagonal density matrix as, This can be obtained by mapping |m → |m m| in Eq. (S9). The purity of this density matrix is equal to γ cl . This classical purity can be used to define a (classical) Renyi-2 entropy for the probability distribution as S cl 2 = − ln γ cl .

Coherence and difference between entropies of the classical and of the quantum models
In a (1 + 1)D classical model, P t (m) represents the probability distribution of 1D states of definite density (labelled by m) of the t th row at time t [26]. When extending such a PCA to a quantum model via (1 + 1)D QCA, the 1D reduced density matrix of the quantum system, ρ t , features, as the diagonal terms, the entries of ρ cl t . The main difference between ρ t and ρ cl t lies in the presence of coherence terms in the density matrix ρ t . As a consequence, the purity of the states ρ t and ρ cl t differ by a contribution which is equal to the 2 -norm of the coherence, C 2 = γ − γ cl , where we have defined γ = Tr ρ 2 t . In terms of C 2 , the difference between the quantum Renyi-2 entropy, S 2 = − ln(γ), and the classical one S cl 2 is = − ln 1 + e S cl 2 C 2 .
(S15) (S16) At criticality, S cl 2 → 0 as t → ∞ and the absorbing state (which is pure) is approached. We can then expand e S cl 2 C 2 as, Keeping only the leading order term and using that C 2 1 when t → ∞ near criticality, further expanding ∆ gives, at first order in the coherence, Therefore, we can approximate, for t 1 and when the system is close to criticality. If, for a given initial condition, the coherence term C 2 (t) decays and approaches zero, we can expect S 2 → S cl 2 . In this case, critical exponents of the quantum and of the classical entropies would agree.

States and Observables
To simulate the evolution of ρ t with MPSs, we use a doubled-space representation of the reduced density matrix. In this framework, operators are mapped to vectors according to the isomorphism |m n| → |m ⊗ |n [39]. One thus has, Here, the symbol ⊗ indicates a product between the two parts of the doubled-space, rather than between rows of the QCA as in the main text and in the previous sections of this supplemental material. Within such a representation, expectation values of observables can be calculated as whereÔ L =Ô ⊗ 1 and |1 = m |m ⊗ |m is the doubled-space vector representation of the identity operator. The purity of ρ(t) can then be calculated as We recall here that the vector representation of the state ρ(t) is normalized in such a way that 1|ρ(t) = 1, so that the purity is in general different from 1. From γ, the second-order Renyi entropy can be calculated as S 2 = − ln (γ).

Time-evolution
In the doubled-space representation, the evolution equation (S2) of the reduced 1D QCA can be expressed as the action of a linear map, Λ t , onto the state |ρ t−1 , i.e., In what follows, the map Λ t is expressed as a matrix product operator, which allows for the application of matrix product state methods to determine the evolution of the matrix product state representation of |ρ t .
To find the form of the map Λ t we can proceed as follows. We first extend the vector representation of |ρ t−1 in order to include also the row t and its doubled-space component. These are both initialized with all particles in the down state. We thus have where the subscripts t − 1, t indicate to which row the vectors entering in the tensor product belong. The second and the fourth entries of the tensor product in the above equation, form the doubled-space components needed to vectorize the density matrix ρ t−1 ⊗ |Ω t Ω t |. The update is obtained, in the density matrix formalism, by applying the global update gate G on both sides of the density matrix and tracing out the degrees of freedom on row t − 1 [c.f. Eq. (S2)]. In the doubled-space representation this is achieved via Here, |1 t−1 is the doubled-space representation of the identity operator with support solely on row t − 1, and implements the partial trace over the associated space. The notation (·) * indicates (element-wise) complex conjugation. Furthermore, we have defined G ij to be the global update which takes, as control sites, the ones described by the vector in the ith entry of the tensor product of Eq. (S23) and, as target sites, those described by the vector in the jth entry of the tensor product. This result shows that the map Λ t appearing in Eq. (S22), acting directly on the vector |ρ t−1 , is given by the operator

TENSOR NETWORK REPRESENTATION OF GLOBAL UPDATE IN DOUBLED-SPACE
In this section we present a method for representing the global update operator Λ t [c.f. Eq. (S24)] as a matrix product operator (MPO). When further representing the state |ρ t as a matrix product state (MPS), standard methods for MPSs can be applied to approximate the time-evolution |ρ t = Λ t |ρ t−1 . These can be found in, e.g., [32][33][34][35].
While this approach is limited to studying reduced dynamics on finite size lattices, it has the significant advantage that it can be integrated into existing tensor network algorithms for dynamics based on MPOs and MPSs, which are extremely common and have been highly optimised.

MPO Representation of G
Before turning to Λ, we construct an MPO representation for G, via the procedure illustrated in Fig. S1, which uses the common diagrammatic notation available for tensor networks [32]. We will consider here the case of three-site neighbourhoods, but this can be generalised. First, G -which is a four-body operator here acting on three control sites and a single target site -is represented as a three-site MPO, see Fig. S1(a). To construct G from these, one simply chooses a particular ordering of target sites and contracts the MPO representations of G for each target site together, according to this ordering, see Fig. S1(b).

MPO Representation of Λ
To evolve ρ t represented in the doubled-space via MPS, we need to find a representation of Λ t [c.f. (S24)] as an MPO in this space. With an MPO representation of G, this can be achieved by repeating the same general procedure for each tensor individually, as illustrated in Fig. S2.
Labelling the i th tensor in a given MPO representation as T [i] , we begin by taking the i th tensor of the MPO representation of G, which we then denote as T i G , and factorising the physical indices, see Fig. S2(a). The tensors representing G |Ω can then be obtained from this by applying a vector representing |• to the appropriate leg, see Fig. S2(b).
The tensors of Λ can then be obtained in two stages. First, a copy of the current tensor is made and the (elementwise) complex conjugate is taken. This new tensor appears in the MPO representation of G * |Ω . This is combined with FIG. S1. Construction of Global Update G as an MPO (a) G is represented as a three-site MPO. For a three-site neighbourhood G is initially represented by a tensor of order-8. Indices relating to the control sites are labelled as l, c, r (left/centre/right) while the target is labelled t. This is then reshaped into a tensor of order-6 by collecting the indices for the central control site and the target site as indicated in the figure. The resulting tensor can be decomposed via, e.g., singular value decomposition [32], to form a three-site MPO. To ensure all local dimensions are equal (corresponding to the vertical legs in the diagram) here we attach the identity matrix, 1, onto the sites to the left and right of the target site in the same row (indicated by t− and t+ respectively). (b) The MPO representation of G is then formed by applying G across the system, once per target site. Here, for illustration, we have chosen a system of size L = 10. The target site for each application of G is highlighted in red. Note that the G obtained has support over 12 sites, which includes a left and a right boundary site. In the main text, these are always assumed to be in the empty state. Since in general G acting on different sites do not commute, the order in which they are applied to form G matters. Here, we have chosen a simple partitioned ordering that minimises the number of layers as shown, though in principle any can be chosen.
the previous tensor by contracting over the legs shown in Fig. S2(c), which correspond to those traced-out via the application of |1 t−1 in the doubled space. The indices of the resulting tensor (which is of order-8) are then collected together as shown to form an order-4 tensor. This is the i th tensor in the MPO representation of Λ t . Repeating this procedure for all i = 1, 2, ..., L + 2 (inclusive of the two boundaries) then produces the desired MPO representation.
To evolve a state |ρ t−1 defined as an MPS over L sites with this MPO in the main text, we first expand |ρ t−1 with an empty site to the left and right resulting in a state with support on L + 2 sites. The MPO for Λ is then applied variationally (see e.g. [34]) before tracing out the left-most and right-most sites to form an approximation of |ρ t .

ERRORS IN ESTIMATION OF CRITICAL EXPONENTS
In this section we discuss the errors for the estimates of the critical exponents of the models. These values are collected in Table S1. For each exponent, the error presented in the main text is taken as the largest of the estimated errors from the various sources that we describe below.

QDKCA model
In determining the critical exponents, one needs to consider errors associated to two primary (though not independent) sources. The first are those associated with estimating the location of the critical point. Such an error is difficult to quantify and can, potentially, be very important. In this regard, however, the QDKCA model holds a key advantage over the QBBR, given that the associated classical DKCA has been studied extensively and the critical point determined to a very high precision. As such, although it remains difficult to quantify the relevance of this error, we can expect this to be negligible, for the QDKCA, as compared to other sources. This is the assumption that we make, and we completely neglect this source of error for the QDKCA.
The remaining error sources for the QDKCA concern the accurate approximation of S 2 and C 2 , and any associated effective exponents. For fixed L, p 1 , p 2 and up to a given t, this error is controlled by the bond-dimension of the MPS, χ. As χ → ∞ this error must vanish. We observe that the tendency of finite-χ errors are to produce an artificial curvature upwards in S 2 (t). This can be seen most clearly in the data of Fig. 2(c) for p 1 = 0.65 between t = [100,200] where, rather than tending to a stationary value, the curve bends upwards. As this curve is not used in approximating the critical exponents, this does not directly affect estimates. However, to avoid under-estimation of errors, the tendency of these finite-χ effects does alter our error analysis indirectly for the QBBR model where finite-χ effects are larger and the critical point is not precisely known. See below for a discussion of this issue.
The value of χ required for an acceptable error can depend strongly on the model's parameters. To estimate the error due to finite-χ, we recalculate each quantity of interest -in this case the critical exponents -using approximations with a χ being half of the one used for the estimate.
The second source of error that concerns the accurate approximation of S 2 and C 2 -which we would like to obtain free of finite-size effects -is that of finite-L. To estimate the relevance of this, we follow a similar procedure as for χ and recalculate quantities using simulations with L/2, but otherwise fixed parameters.
Finally, an additional issue when estimating the values of critical exponents are finite-time errors. To quantify these, we recalculate the estimated quantities (via fits or averages) but over an interval of time that is half of that used for the original estimate, starting at the same initial time.
As discussed in the main text, q ent is estimated via a power-law fit of the form t −qent to S 2 (t) approximated with χ = 128, L = 256 between t = [50, 200] for p 1 = 0.645. This gives q ent = 0.21265 (all estimated quantities are given in this section to five decimal places while estimated errors are given to three significant digits). Repeating this procedure with S 2 (t) calculated using χ = 128, L = 128 gives q ent = 0.21270. With χ = 64, L = 256 we obtain q ent = 0.21263. Finally, with χ = 128, L = 256 between t = [50, 125] gives q ent = 0.21635. Taking the absolute difference between each of these estimates and the original, we approximate the errors due to each source to be: 5.37 × 10 −5 , 2.21 × 10 −5 , and 3.70 × 10 −3 for finite-L, finite-χ and finite-time respectively.
To estimate errors in q coh the same procedure is used as with q ent for error due to finite-χ and finite-time. In this case, the original estimate is obtained by averaging over the effective exponent q coh (t) from t = [450, 550] to obtain q coh = 2.96486. Repeating the error calculations as for q ent obtains: 2.88 × 10 −3 for finite-χ errors and 3.57 × 10 −3 for finite-time errors (this being obtained by averaging over t = [450, 500]). Due to the long times required, estimating errors due to finite-L from calculations with L = 128 leads to a significant overestimation, as for this time the state is already approaching the absorbing state. As such, to estimate finite-L errors in this case we instead compare the results with those obtained from a method, free of finite-size effects, which provides an estimated error of 1.3 × 10 −3 . Details of this method can be found in [40]. These data were produced using χ = 64 but otherwise equal model parameters.

QBBR model
While the error due to the determination of the critical point was assumed to be minimal for the QDKCA, in the QBBR model -since the classical BBR model is much less studied -this can play an important role. This is particularly true since the estimation of critical exponents is very sensitive to the distance from the true critical point.
For APTs, a standard procedure for estimating errors in this regard is to bound the value of a critical exponent using curves that are known to be in a given phase [4]. In the case of a power-law decay at a critical point, curves in the inactive phase (which decay more rapidly than the critical curve) can provide upper-bounds to the exponent. Similarly, curves in the active phase can provide lower-bounds. By simulating a set of curves on some grid up to a given time, one then identifies curves in a given phase and uses these for the error estimate. To decrease the error in such an estimate, a finer grid can be used. However, distinguishing the phase of a given curve will require longer and longer times the closer it is to the critical point. As such, ultimately, the limitation in this procedure stems from the maximum time that can be accurately approximated.
For the QBBR model, we estimate errors due to the uncertainty on the critical point using a simple grid search. That is, from the values of p 1 chosen -which form the "grid" -we select as the estimated critical point the one that results in a curve S 2 (t) that best approximates a power-law. Errors are then estimated by recalculating quantities using two other curves, one from the inactive phase, and one from the active phase. For the curves of S 2 (t) presented in the main text [c.f. Fig. 2(c)] the one with p 1 = 0.61 best approximates a power-law and is therefore chosen as the estimated critical point. The curve below this one, with p 1 = 0.6, shows clear behaviour characteristic of the absorbing (inactive) phase. As such, we take this curve to provide an upper bound on q ent . To provide a lower-bound, the curve with p 1 = 0.62 displays curvature indicative of belonging to the active phase. However, as mentioned previously, finite-χ errors tend to produce this an upward curvature artificially. As such, to avoid underestimation of this error, we instead take the curve with p 1 = 0.625 in this instance.
The estimates of q coh for the QBBR model proceed in a similar fashion. Taking the same p 1 = 0.61 to estimate the critical point and averaging over t = [150,200] gives q coh = 2.73746. The finite-L, finite-χ and finite-time errors (calculated by averaging over t = [150, 176]) are estimated as 7.11 × 10 −3 , 1.56 × 10 −3 and 3.75 × 10 −3 respectively.
For estimating the errors related to determination of the critical point, only values of p 1 for which q coh (t) is approximately constant over a substantial interval of t can be used. In practice, we find that this excludes values of p 1 from the active phase as the effective exponent tends to diverge, meaning no such interval can be found. As such, here we construct a simple estimate using the curve that previously provided the lower bound for q ent (i.e. the curve with p 1 = 0.60) and assume a symmetric error about the estimate of q coh . Calculating q coh just as for p 1 = 0.61 but with p 1 = 0.60 provides q coh = 2.47247, giving an associated error of 2.65 × 10 −1 .
The error estimates from the various sources discussed are summarised in Table S1.