Rule 54: Exactly solvable model of nonequilibrium statistical mechanics

We review recent results on an exactly solvable model of nonequilibrium statistical mechanics, specifically the classical Rule 54 reversible cellular automaton and some of its quantum extensions. We discuss the exact microscopic description of nonequilibrium dynamics as well as the equilibrium and nonequilibrium stationary states. This allows us to obtain a rigorous handle on the corresponding emergent hydrodynamic description, which is treated as well. Specifically, we focus on two different paradigms of Rule 54 dynamics. Firstly, we consider a finite chain driven by stochastic boundaries, where we provide exact matrix product descriptions of the nonequilibrium steady state, most relevant decay modes, as well as the eigenvector of the tilted Markov chain yielding exact large deviations for a broad class of local and extensive observables. Secondly, we treat the explicit dynamics of macro-states on an infinite lattice and discuss exact closed form results for dynamical structure factor, multi-time-correlation functions and inhomogeneous quenches. Remarkably, these results prove that the model, despite its simplicity, behaves like a regular fluid with coexistence of ballistic (sound) and diffusive (heat) transport. Finally, we briefly discuss quantum interpretation of Rule 54 dynamics and explicit results on dynamical spreading of local operators and operator entanglement.


Introduction
Exactly solved models are a major cornerstone of statistical mechanics and physics in general. While free (quadratic) models and their perturbations provide some insight, to achieve realistic statistical physical behaviour found in real materials requires strongly interacting models of which only several solvable classes are known. The first and most notable class, in the context of equilibrium physics, are two-dimensional lattice models related to solutions of the Yang-Baxter equation which give crucial exact information about the universality classes of statistical systems in two dimensions [1], and relate to Bethe-ansatz solvable quantum models in one dimension [2][3][4]. However, for studying out of equilibrium properties, in particular for time-dependent correlation functions or quantum quenches, such Yang-Baxter-Bethe solvable models represent a much harder challenge. Computation of time correlation functions through the so-called formfactor-expansion is a formidable task, so far accomplished with only partial success in particular models [5][6][7][8][9][10][11]. Nevertheless, integrability techniques resulted in a successful hydrodynamic approach in such models. The most notable of which is the so-called generalised hydrodynamics (GHD) [12][13][14], which has achieved remarkable success in predicting large space-time scale behaviour of observables in integrable systems.
However, hydrodynamics is not a rigorous theory and in particular it relies on the assumption of a clear separation of space-time scales. The fact that sub-ballistic corrections generically result in diffusion, as derived within GHD in Ref. [15] (see also [16]), could be at least partly attributed to the central assumption behind the hydrodynamic picture, which is the immediate loss of memory (correlations) in the quasi-particle scattering processes. It is thus of utmost importance to have at our disposal another type of model -or a class of models -with generic physical behaviour and for which dynamical physical quantities are accessible by a rigorous analysis free from assumptions. Within such a class of models we can then achieve the 'holy grail' of nonequilibrium statistical physics, which is to derive macroscopic transport laws from reversible and deterministic microscopic equations of motion [17].
Indeed, in the last several years it has been recognised that such a model exists and seemingly belongs to a class distinct from regular Bethe-ansatz solvable models. It is the Rule 54 of the family of reversible cellular automata (RCA54) proposed and classified in Ref. [18]. The model RCA54 can also be understood as a staggered sublattice version of the model 250R of reversible local automata as classified earlier by Takesue [19], or a discrete-time deterministic limit of the Fredrickson-Andersen model [20]. It is arguably the simplest microscopic (deterministic) physical theory in 1+1 dimensions with strong local interactions and asymptotically freely propagating excitations -quasiparticles, or solitons. RCA54 has been proposed to be integrable already in [18] based on mainly qualitative arguments. The first exact solution, however, came in Ref. [21], where the nonequilibrium steady state of the model driven by a pair of chemical baths at the boundaries was analytically found. This led to many other results, such as generalisation to larger families of boundary driving [22,23], diagonalisation of the Markov propagator [23,24], and exact large-deviation results in the boundary driven setup [25]. Later, many properties of the model on an infinite lattice without stochastic boundaries have been exactly obtained, such as dynamical structure factor and multitime correlation functions [26][27][28]. The model has also been studied in the quantum context [29], and its simple but non-trivial dynamics was found to provide an ideal setting to study large-scale properties of the physics of local observables and operator spreading [30][31][32][33][34][35][36][37]. Furthermore, an intriguing connection with the dynamics of T Tdeformed conformal field theories was recently established [38,39].
In particular, with respect to transport of quasiparticle excitations (or solitons), Rule 54 provides a model of a generic physical fluid with coexistence of ballistic (convective) and diffusive (conductive) transport. Remarkably, unlike typical Betheansatz solvable systems, this model in many instances allows for fully closed-form solutions despite being interacting and thus it allows the understanding of the aforementioned generic transport behaviour on a microscopic level. The purpose of this review is to provide a comprehensive overview of recent results on the Rule 54 model and discuss their comparison with simple predictions of hydrodynamic theory.

Definition of the model and summary of the results covered
We consider deterministic dynamics defined on a 1 + 1 dimensional discrete lattice with space-time points labelled by (x, t) ∈ Z 2 and a field variable s t x taking only binary values s t x ∈ Z 2 = {0, 1}. We may restrict the dynamics only to a staggered (diamond, or light-cone) sublattice of Z 2 of points (x, t) satisfying x + t = 0 (mod 2) and define s t+1 for some binary function χ : Z 3 2 → Z 2 . Specifically, the function with the binary code 54 reads s 2 = χ(s 1 , s 2 , s 3 ) ≡ s 1 + s 2 + s 3 + s 1 s 3 (mod 2), and can be graphically represented as where white/black boxes represent field values 0/1, and red-framed boxes correspond to the updated values at the next time instance. As the graphical representation suggests, the middle site is changed whenever at least one of the neighbours is black. Dynamics is generated by two sequences of parallel updates (see Fig. 1), which in two steps maps a configuration over a zig-zag saw s t−1 first to s t and then to s t+1 , where s t ≡ (. . . , s t−1 x−1 , s t x , s t−1 x+1 , . . .). An example of a typical trajectory can be seen in Fig. 2 (ignoring the boundaries for the time being). The dynamics can be interpreted as a gas of solitons travelling with speeds ±1 which scatter pair-wise, while each scattering produces a time-lag (or shift) of a soliton of one step in time (or space). In particular, the diagrams in (3) can be interpreted as different possible local configurations of solitonic dynamics. The first diagram represents empty space, diagrams 2 and 7 show a left mover, diagrams 4 and 5 represent a right mover, while the rest of the diagrams correspond to different steps of the scattering event between two oppositely-moving solitons-the two solitons first merge, temporarily disappear and reappear immediately afterwards (given by diagrams 6, 3 and 8 respectively).
The key problem discussed in this review is how to achieve the full understanding of equilibrium and nonequilibrium statistical physics of this model. This goal can be pursued within two paradigms of statistical mechanics. In the first, we consider the system defined on a finite lattice Z n , of even size n, and we couple the left and the right edge of the chain to stochastic baths of solitons. After removing the degrees of freedom of the reservoirs we end up with a perfect Markov chain model, where the cells in the bulk are updated deterministically, while the cells near the boundaries are evolved stochastically. The set of parameters characterizing the boundary driving uniquely determines the nonequilibrium stationary state (NESS) that the system approaches at long times. In Section 3 we show how to provide exactly solvable Markovian boundaries and how to construct a correlated NESS in terms of the so-called patch-state ansatz. Additionally, we prove a general ergodicity theorem for the boundary driven setup, which guarantees the uniqueness and exponential relaxation to the steady state. In Section 4 we then provide a solution to the steady state cancellation mechanism in terms of a simple cubic algebra whose representation yields a compact matrix product solution to the full steady state. Remarkably, this form can be extended to subleading eigenvectors of the Markov propagator and allows us to obtain an essential part of its spectrum, which, in particular, gives access to the spectral gap characterising the relaxation to the NESS. Moreover, the same cubic algebra can be further generalised to generate an exact large deviation function describing fluctuations (in the steady state) of a large class of spatially inhomogeneous observables, as elaborated in Section 5.
The second main paradigm considers the dynamics of the thermodynamic states of the system defined on an infinite lattice, or equivalently, the dynamics of local observables of the large but finite system far from the boundaries. In Section 6 we provide a basic hydrodynamic description of the model using the two elementary local conserved charges: the densities of left and right movers. In Section 7 we provide a space-time dual description of the model. We start by considering an efficient matrixproduct-state representation of stationary probability distributions of configurations in time, which determine all multi-time correlations of ultra-local observables, evaluated in equilibrium thermodynamic states. We then proceed by showing that RCA54 allows for a deterministic description also when considering the space-evolution of timeconfigurations. For facilitating computations we introduce efficient diagrammatics, akin to quantum circuit diagrams, which are used to compactly encode a diversity of rather formal algebraic relations. Arguably the strongest result, elaborated in Section 8, is the construction of an exact time-dependent matrix product ansatz (tMPA) for time evolution of arbitrary local observables. Although the dimension of tMPA is formally infinite, it effectively grows only as ∝ t 2 when time evolution up to time t is considered. The result allows for a number of explicit computations, such as the dynamical structure factor (two-point space-time correlation function) and the timedependent density profiles following an inhomogeneous quench. Furthermore, it provides a formal bound on the rate of operator spreading for the quantum interpretation of the model.

Boundary driven cellular automaton and general equilibrium states
We begin by considering the model subject to stochastic boundary driving. The system is defined on an even number of sites n with instantaneous dynamical variables specified over a zig-zag saw , s = (s 1 , s 2 , s 3 . . . , s n ) ≡ (s 1 1 , s 0 2 , s 1 3 . . . , s 0 n ) ∈ C n , where C n = Z ×n 2 is the set of all configurations on the lattice of size n. The macroscopic state of the system p is a probability distribution over the set of all configurations and can be interpreted as a 2 n -dimensional vector with p s ≥ 0 denoting the probability of configuration s. If dynamics were completely deterministic, such a probabilistic description of a finite system would be redundant. However, we allow for stochastic updates of boundary pairs of cells (s 1 , s 2 ), and (s n−1 , s n ) which cannot be consistently fixed by the deterministic rule (2). We make a minimal assumption that the stochastic updates of boundary cells are local and Markovian, i.e. the probability of jumps of s 1 (or s n ) can only depend on cells s 1 , s 2 (or s n−1 , s n ) in the immediate space-time neighbourhood. We thus define the time evolution of the state vector p(t) starting from some initial state p(0) as, x t Figure 2. A snapshot of Monte Carlo dynamics of nonequilibrium stochastically boundary driven deterministic CA, Rule 54, for n = 80, α = 0.9, β = 0.1, γ = 0.4, δ = 0.6 (Bernoulli driving (13) with ζ = η = 1/2). Time runs downwards. Grey squares denote occupied environmental cells which are generated by a Bernoulli shift with probability 1/2, while blue (red) squares are occupied boundary cells determined via ultralocal Markov chains. Note that the right end is "hotter" than the left one and that the average (steady-state) soliton current flows to the left.
where the one-period propagator U (2 n × 2 n matrix) is factored in two half-time steps, as shown in Fig. 3, Each step is now given in terms of parallel updates of even or odd sites where U j−1 j j+1 = 1 2 j−2 ⊗U ⊗1 2 n−j−1 , with U being an 8×8 matrix, acts non-trivially on a triple of adjacent sites (j − 1, j, j + 1) and only affects the cell at position j depending on the values of cells at positions (j − 1, j, j + 1) according to the Rule 54. Specifically: Figure 3. Illustration of the stochastic boundary driving (7) with two site boundary stochastic maps U L,R (see Eqs. (8,9)). In blue/red we denote the sites before/after the update.
The boundary updates are generated in terms of 4 × 4 stochastic matrices acting on the left-most (right-most) two sites. As U L,R are stochastic, i.e. their non-negative matrix elements in each column add to one, the full propagator U is a stochastic matrix conserving probability (see Fig. 3). The boundary matrices have to satisfy the following compatibility conditions essentially implying that only the boundary cell may get stochastically updated, while the site next to it (i.e. n − 1 in U R n−1 n and 2 in U L 12 ) acts as a control cell, analogously to the non-central cells in the bulk propagator.
Requiring that the many-body Markov chain process (6) admits an exactly solvable nonequilibrium steady state (NESS) -an eigenvector p 0 of eigenvalue 1, Up 0 = p 0puts additional constraints to the boundary operators U L , U R . A complete classification of exactly solvable local (2-site) boundaries is composed of two multi-parametric families. The first is the so-called Bernoulli driving (originally proposed in [21] and generalised in [22]): where Each bath is parametrised by a triple of independent parameters: ζ, α, β ∈ [0, 1] for the left boundary, and η, γ, δ ∈ [0, 1] for the right boundary. They admit a simple interpretation in terms of a composition of an ultralocal Markov chain (a 2×2 stochastic matrix parametrised by probabilities α, β on the left and γ, δ on the right), and a local Rule 54 map centred around the same site, while the value of the imaginary cell inside the bath is given in terms of a Bernoulli process (with the probability ζ on the left and η on the right). Another distinct family is the conditional driving [23]: where α, β, γ, δ ∈ [0, 1] are some driving rates parametrising the left and the right bath. We refer to this as conditional driving, since in U L 12 (U R n−1 n ) the probability of changing the site 1 depends only on the state of the neighbouring site 2 (changing the site n depends only on the state of the site n − 1). For instance, if the site 2 is in state 0 then the site 1 will be stochastically set to state 0 with the rate α or to state 1 with the rate 1 − α. On the other hand, if the site 2 is in the state 1, the site 1 will be set to state 0 or 1 with the rates β or 1 − β, respectively. The analogous holds for U R n−1 n . These two classes of boundary driving maintain the integrability of the model, which may be intuitively understood as a consequence of the fact that they (with some probabilities) create or destroy solitons at the boundaries. We note that the information about three-site subconfigurations is needed to completely characterize the local particle content, therefore a minimal matrix product ansatz providing such a state should admit a 3-dimensional auxiliary space. This fact will be important in the sections below, when dealing with solutions to boundary driving.
Proof. According to the Perron-Frobenius theorem [40], a finite, non-negative matrix U is irreducible, if for any pair of configurations s, s ∈ C n there exists a natural number t 0 ∈ N such that the corresponding matrix element of the t 0 -th power is nonzero, An irreducible matrix U is aperiodic if for any configuration s ∈ C n , the greatest common divisor of the set of recurrence times {t j } is 1,  Illustration of the proof of irreducibility and aperiodicity of the Markov matrix U. Blue and red configurations, s = (1, 0, 0, 1, 1, 1, 1, 0, 1, 0) and s = (0, 0, 1, 0, 1, 1, 0, 0, 0, 1), are connected with the Markov-graph walk for generic probabilities 0 < α, β, γ, δ < 1 in at least t 0 = 15 time steps where the boundary cells are chosen as indicated by green cells (the boundary conditions are generated by causal/anti-causal absorbing boundaries in the upper/lower part of the walk). The values of the boundary cells are thus determined by copying the values of the near-by bulk cells in the direction of the grey arrows. Consequently [U t0+tgap ] s,s > 0 for any t gap ≥ 0 (and any other pair of initial/final configurations s, s with possibly different t 0 ), which implies irreducibility and aperiodicity of U.
Let us first show irreducibility. By the definition of boundary propagators (cf. (13) and (15)), the Markov matrix U connects each configuration s to exactly 4 other configurations s , which exhibit all the possible configurations of boundary bits, (s 1 , s n ) ∈ {(0, 0), (0, 1), (1, 0), (1, 1)}. This holds for all values of parameters α, β, γ, δ (and ζ, η), except for the marginal case when some of the parameters are equal to 0 or 1, but this option is excluded by the assumption of the theorem. Let us now take a sufficiently large positive integer t 0 , to be determined below, and fix s(t 0 ) ≡ s , s(0) ≡ s. We shall then construct a walk which can be understood as a path through the Markov graph defined by positive elements of U that connects s and s in t 0 steps and implies [U t 0 ] s ,s > 0 (see Fig. 4 for a 'self-contained' graphic illustration of the idea of proof). We are still free to choose the values of the boundary cells s 1,n (t) along the walk t ∈ {1, 2, . . . , t 0 − 1} apart from the ends. For the first part of the walk t = 1, 2 . . . t + , up to some t + < t 0 , we are fixing them with the rule The evolution of the interior values of the cells s x (t) for 1 < x < n and t ≤ t + is then completely specified by the deterministic RCA54, while (19) provide the causal absorbing boundary conditions. Indeed, each time the boundary cell, say x = 1, gets occupied, s 1 (t) = 1, the soliton is absorbed (see Fig. 4). As the solitons only move ballistically (at speed 1) and scatter pairwise (with time-lag 1), it is clear that a finite time scale t + ∈ N exists, surely smaller than n 2 , after which all the solitons will be absorbed and we end up in a vacuum configuration s(t + ) = (0, 0 . . . , 0). For the rest of the walk t ∈ {t + +1, . . . t 0 } we need to show that alternative boundary rules exists, which create the configuration s out of the vacuum in another t − = t 0 − t + steps. This is easily achieved by using time-reversibility of RCA54 and arguing that a vacuum configuration is again generated from s in some t − steps if the anti-causal absorbing boundary conditions are set, which are equivalent to (19) when the time runs backwards, The entire walk then connects s to s in t 0 = t + + t − steps and implies [U t 0 ] s ,s > 0, for an arbitrary pair s, s ∈ C n , with the minimal possible integer t 0 in general depending on the choice of s, s . This proves irreducibility of (7).
Considering s = s, we have just shown that [U t 0 ] s,s > 0 for some t 0 = t + + t − depending on s. However, in between annihilating the configuration s in t + time steps and then creating it again in another t − steps ‡, we can stay in the vacuum state for an arbitrary additional number of steps t gap ≥ 0. In this way the walk is increased by a segment of t gap intermediate vacuum configurations, and we still have a non-zero matrix element, The greatest common divisor of the set {t 0 + t gap ; t gap ∈ Z + } is clearly 1, so we have shown aperiodicity.
In conclusion, the Perron-Frobenius theorem [40] guarantees that the NESS probability state vector p 0 , satisfying the fixed point condition Up 0 = p 0 , is unique and all other eigenvalues of U lie strictly inside the unit circle. As a consequence, the Markov dynamics (6) is ergodic and mixing and an arbitrary initial probability state vector p(0) converges to NESS exponentially fast in t.
p' Figure 5. Illustration of the patch state ansatz (23) for NESS probability state vectors p, p .

NESS: Patch State Ansatz
The NESS probability state vector p = p 0 can be split into even and odd time slice, satisfying a pair of fixed point equations We shall now explicitly construct the probability state vectors p and p , by solving Eq. (22) in terms of a simple ansatz, which we term a patch state ansatz (PSA) [illustrated in Fig. 5].
For either of the two boundary driving families (13,15), the NESS solution p, p ∈ R 2 n of the fixed point condition (22) can be written in the form p s 1 ,s 2 ,...,sn = L s 1 s 2 s 3 X s 2 s 3 s 4 s 5 X s 4 s 5 s 6 s 7 · · · X s n−4 s n−3 s n−2 s n−1 R s n−2 s n−1 sn , p s 1 ,s 2 ,...,sn = L s 1 s 2 s 3 X s 2 s 3 s 4 s 5 X s 4 s 5 s 6 s 7 · · · X s n−4 s n−3 s n−2 s n−1 R s n−2 s n−1 sn , for some rank-4 and rank-3 tensors of strictly positive components X ss uu , X ss uu , L suu , L suu , R ss u , R ss u , with binary indices s, s , u, u ∈ {0, 1}. To uniquely determine the algebraic expressions for these tensors in terms of the parameters of the model, one has to plug the ansatz (23) into the NESS condition (22). First we find a minimal sufficient set of equations that the tensors have to satisfy, by fixing the normalization and taking into account the gauge symmetry of the ansatz. The normalization of the PSA (23) can be chosen such that X 0000 = X 0000 = 1, L 000 = R 000 = 1.
The total number of 2 × 16 + 4 × 8 − 4 = 60 unknowns can be further reduced by exploiting the following gauge symmetry which conserves the patch ansatz (23), as well as the defining equations (25,26) for arbitrary nonzero gauge 'fields' f ss , g ss .
While a detailed analysis and explicit expressions for the PSA tensors in terms of boundary parameters can be found in references [21,22], here we only discuss a generic solution form of the bulk part of the equations. Specifically, the bulk equations (25) can be solved independently (before incorporating boundary conditions), and the solution can be parametrised uniquely -up to a choice of gauge (27) -in terms of two free (spectral) parameters ξ, ω X ss tt = T (ss ),(tt ) (ξ, ω), X ss tt = T (ss ),(tt ) (ω, ξ), T (ξ, ω) = The boundary equations (26) then fix the spectral parameters ξ, ω as functions of the boundary driving parameters α, β, γ, δ (and ξ, η) [see Section 4]. The remarkable fact is that the solutions do not explicitly depend on the system size n, hence it is clear that such a NESS can only describe ballistic transport (i.e. net soliton current independent of n).

Conserved charges
The 4 × 4 matrix T (ξ, ω), defined in Eq. (28), can be interpreted as a two-parametric transfer matrix generating the local charges of the RCA54 model. To simplify the discussion, let us assume periodic boundary conditions, n + 1 ≡ 1, and consider a steady state vector p(ξ, ω) whose components are given in terms of the transfer matrix as, For any pair of positive real parameters ξ, ω > 0, p(ξ, ω) represents the statistical ensemble -state that is invariant under time translation (as well as translationally invariant in space) whose partition sum is given simply in terms of powers of the transfer matrix as Parameters log ξ, log ω can be interpreted as chemical potentials corresponding to the conserved numbers of left and right movers, N − and N + , Specifically, the latter can be generated directly from the logarithmic derivatives of the transfer matrix which can clearly be written as extensive sums of local densities These are the elementary local charges of RCA54 and will be discussed later in the context of hydrodynamics of the model. We note that N ± by no means represent a complete set of local charges. Specifically, searching for all local charges with densities of support size > 3, we find of the order of Fibonacci number F = F −1 + F −2 of trivially conserved charges corresponding to F non-interacting unidirectional soliton configurations. See the discussion in [41] for the details.

Matrix product ansatz and Markovian excitations
In this section we will recast and generalise the exact solution for the NESS in terms of a more standard matrix product ansatz (MPA). From now on we fix the following convention: notation for row/column vectors in the auxiliary space V a will be Dirac bras/kets. Quantities that are vectors in physical space will be denoted by bold-face Roman letters as before, e.g. the probability state vector p(t). The numeral subscript of an operator or vector in the physical space is the site position in the tensor product physical space (R 2 ) ⊗n on which it acts nontrivially. The physical space components will be labelled by a binary index and written in corresponding non-bold font, e.g. p s 1 ,...,sn . Matrices (operators that act nontrivially either in auxiliary or physical space) will be denoted by capital Roman letters. We are interested in solving the eigenvalue problem for the Markov propagator, which can be conveniently split into two coupled linear equations for the even and odd half-steps of the period of time propagation (generalising (22)), Once we solve this, we have access to the full dynamics of the probabilities at time t in terms of eigenvalues Λ j (Λ 0 = 1) and the corresponding eigenvectors p j (assuming that U is not defective), where c j are constants given by the initial probability distribution p(0). Using this labelling, the eigenvector p 0 denotes the NESS as it does not decay in time, while the components along the rest of the eigenvectors p j , j ≥ 1, decay exponentially with the rate − log |Λ j |. If all multiplicative rates are bounded away from the unit circle |Λ j | < 1, any initial state p(0) will relax to the NESS. In principle, one may find eigenvalues that lie on the unit circle, but are not 1. These would correspond to persistently oscillating eigenvectors analogous to the ones for quantum Markov processes [42][43][44][45]. In our case, however, their existence is prohibited by the Perron-Frobenius theorem (see the discussion in Subsection 3.1).

NESS: Matrix Product Ansatz
Following [23,25] we take a staggered matrix product ansatz (MPA) for the NESS in the form, We will use this as an ansatz to solve the split eigenvalue equation (36). We require that the operators W s , W s and some operator S satisfy the following bulk algebraic relation, Component-wise the above vector equation reads A representation of this algebra can be found on a 3-dimensional auxiliary vector space V a = C 3 in terms of matrices while the other set of auxiliary space matrices W s (ξ, ω) is obtained from W s (ξ, ω) by swapping the parameters, The representation is parametrised in terms of two (so far) free variables ξ and ω, which (as we argue later) correspond precisely to the variables parametrising the transfer matrix introduced in (28), therefore we shall refer to them as spectral parameters. We note that S 2 = 1, which together with U 2 = 1 and the bijection between the two sets of matrices implies a dual bulk relation, We now turn to the boundaries. We demand that the spectral parameters ξ, ω, eigenvalue parameters Λ L , Λ R and the boundary vectors l 1 |, l 12 |, |r 12 , |r 1 , satisfy the following set of equations, If these equations are satisfied, the pair of states p, p written in terms of MPAs (38) and (39), solves the staggered eigenvalue equations (36). This can be straightforwardly verified by plugging the MPAs into e.g. the first equation of the set (36), and apply the appropriate relations to transform p into p , U e p = U 123 · · · U n−5 n−4 n−3 U n−3 n−2 n−1 U R n−1 n l 1 |W 2 W 3 · · · W n−4 W n−3 W n−2 |r n−1 n = Λ R U 123 · · · U n−5 n−4 n−3 U n−3 n−2 n−1 l 1 |W 2 W 3 · · · W n−4 W n−3 W n−2 W n−1 S|r n = Λ R U 123 · · · U n−5 n−4 n−3 l 1 |W 2 W 3 · · · W n−4 W n−3 SW n−2 W n−1 |r n = · · · = Λ R U 123 l 1 |W 2 W 3 SW 4 W 5 · · · W n−1 |r n = Λ R l 12 |W 3 W 4 W 5 · · · W n−1 |r n = Λ R p .
To get to the second line, we apply the boundary propagator U R n−1 n , which, according to relation (47), replaces the right boundary vector |r n−1 n with |r n , and at the same time introduces the delimiter matrix S. We continue by applying the dual bulk relation (44), which moves S to the left, while the roles of W s and W s are exchanged (fourth line). We finish by absorbing the matrix S at the left and finally change l 1 | into l 12 |. The second part of (36) is proven analogously.
Solving the left boundary equations (45,46) for ξ, ω and boundary vectors l s | , l ss | (given explicitly in Appendix A) we obtain Similarly, solving the right boundary equations (47,48) for ξ, ω and the |r s , |r ss (again listed in Appendix A) we get These and all explicit result below are written for the specific case of conditional driving (15), while results for the Bernoulli family are obtained fully analogously [23]. Equating ξ and ω from these pairs of equations yields an eigenvalue equation, which expressed in terms of the eigenvalue Λ = Λ L Λ R reads Importantly, Λ = 1 is always a solution. The corresponding eigenvector is the NESS. The remainder is a cubic polynomial giving three other eigenvalues corresponding to three decay modes. These eigenvalues do not depend on the system size n. We refer to these four eigenvalues as the NESS orbital.

Markovian excitations
In order to generalise the results from the previous subsection to other eigenvalues and eigenvectors (decay modes) of the Markov matrix U we will follow the spirit of the standard coordinate Bethe ansatz (see e.g. [46]) for finding solutions to eigenstates containing interacting particles. However, we need to adapt it to a coordinate matrix product ansatz picture similar to the one used for the asymmetric simple exclusion processes [47]. As we will see the leading decay modes can be understood in terms of localized particle excitations of the NESS, the latter being analogous to the vacuum state. The leading decay modes will thus be a superposition of single-particle excitations. The leading decay mode is the eigenvector of U corresponding to the eigenvalue Λ 1 with the real part closest to 1, Apart from the eigenvectors with eigenvalues on the unit circle, it is the most dynamically relevant eigenvector, since it dominates the long-time asymptotic relaxation. In this subsection, we will provide a compact MPA representation of the eigenvectors corresponding to the orbital that contains the leading decay mode. The generalisation to other orbitals remains an open question, although the complete set of eigenvalues can be predicted by a simple conjecture. The starting point is to double the auxiliary space, which now consists of two copies of the original auxiliary space V a = V a ⊕ V a = C 2 ⊗ V a . We generalise the W operator to the doubled auxiliary space V a as, where e ij = |i j|, i, j ∈ {1, 2} is the Weyl basis of 2 × 2 matrices. We define a set of where we introduced the operators F ( ) only on a single copy of the original auxiliary space. The explicit representation of these operators is given in Ref. [23] in terms of a four dimensional auxiliary space V a = C 4 . Note that unlike the MPA solution corresponding to the NESS orbital, we have so far been unable to find a 3-dimensional representation of auxiliary space operators, so a 4-dimensional single-copy auxiliary space is assumed here. The operators (55) linearly depend on two complex coefficients c + , c − , and in the inhomogeneous MPA (introduced below), the complex variable z plays the role of the multiplicative quasi-momentum of the quasiparticle excitation on top of the NESS, which are created by the operators (55).
We now introduce the following inhomogeneous (site-dependent) MPA and p is the eigenvector of U = U e U o (8,9) with eigenvalue Λ = Λ L Λ R , like before. The interpretation of z as the momentum of a single excitation should now be clear. Due to the presence of the off-diagonal nilpotent e 12 in the site dependent matricesF (k) ,Ĝ (k) these can only create at most one excitation on the NESS in the MPA (56) and (57). Furthermore, since the site dependent matrices contain both terms with z and z −1 , this MPA is like a standing wave. We deform the bulk relation and its dual, (40,44), to an inhomogeneous form -where a free integer k denotes a spatial coordinate -which our newly introduced operators have to satisfy whereŜ = 1⊗S. Analogously as before we demand the boundary vectors and the other parameters to provide a nontrivial solution to a set of generalised boundary equations, To understand why p is indeed an eigenvector, let us for clarity consider an example with fixed n = 8 spin sites, (4)Ŵ 6 |r 78 = . . . We start by first acting with U 123 (note that all the local operators commute, cf. (12), and can be applied in any order), which using (60) 3Ŝ to the penultimate right hand side via the bulk relation (58), . . Prior to acting with the final U n−3 n−2 n−1 , we apply U R so it createsF (n−3)Ŵ n−1Ŝ at the right end, . . We then use the top Eq. (58) for the final U n−3,n−2,n−1 to moveK (n−5)Ŵ n−3Ŝ to the far end intoK (n−3)Ŵ n−1Ŝ and we thus finally obtain the form (57) of p , The odd-part of the propagator U o acts analogously in reverse, hence we obtain the second part of the eigenvalue equation, Note that for z = 1, c + = c − = 0, one recovers the NESS bulk relations (40,44). The boundary equations given by (59)(60)(61)(62) can be solved to obtain the parameters ξ, ω, z, c + , c − ∈ C, together with the eigenvectors and eigenvalues of the first (leading decay-mode) orbital. Explicitly, from the left boundary equations we obtain and from the right ones, where m = n 2 − 2. Pairwise identifying expressions for ξ, ω and c + /c − from the left and right boundary equations, we obtain a triple of coupled equations for the unknowns Λ L , Λ R and z. For the details of the solution procedure and the form of the boundary Figure 6. The spectrum of U for n = 12 and α = 1/2, β = 1/3, γ = 1/5, δ = 4/7. The black crosses show numerical results. The red (p = 1), green (p = 2), brown (p = 3), orange (p = 4) points are solutions for the p-th orbital eigenvalues (68). The blue squares are the roots of characteristic polynomial for the NESS orbital. Note that eigenvalues have to be largely degenerate on average as the number of distinct eigenvalues is 4+18×8 = 148 ∝ n 2 , while the total number of eigenvalues is 4096 = 2 n . The dashed grey curve is the unit circle. The blue curve is an algebraic curve to which the leading decay mode orbital converges in the thermodynamic limit.
vectors we refer the reader to [23]. Using this solution it can be shown that the eigenvalue of the leading decay mode Λ 1 scales with large n as The 1/n scaling of the gap implies that the relaxation time is proportional to n, which is consistent with the ballistic transport observed in [21].
Based on these results one can formulate a conjecture for the higher orbital eigenvalues [23,24]. Let p denote the orbital number, 1 ≤ p ≤ m = n/2 − 2. Then the complete set of eigenvalues Λ = Λ L Λ R satisfy the following set of Bethe-like equations, where p = 1 corresponds to the first orbital obtained analytically via Eqs. (65)(66).
The eigenvalue solutions of (68) are shown in Fig. 6 and agree perfectly with numerical diagonalization of the Markov matrix U. In contrast to the usual Bethe ansatz [2], this conjecture implies that we have only one momentum parameter (z) even for manyparticle excitations. This is consistent with the Bethe-ansatz description of Rule 54 proposed in [32], and can intuitively be understood as a consequence of a very singular dispersion relation, namely all quasi-particles (solitons) move with the same speed ±1. Another many-body effect is a large degeneracy of each eigenvalue in the orbital p which is conjectured [23] to be exponential, specifically 2 p−1 .
Expanding the last equation (68) for in the leading order of 1/n, we get that z becomes unimodular, z = e iκ , where the momentum κ is restricted to the interval [0, π). Note that if some eigenvalues Λ L Λ R in the n → ∞ limit converged to points on the unit circle different from 1, this would imply persistent oscillations and/or non-decaying fluctuations in time for the model in the thermodynamic limit. Although one may naively expect that the effect of the boundaries becomes thermodynamically irrelevant, pushing all eigenvalues to collapse onto the unit circle (similarly to the closed model studied in Sections 6-8), this is not the case. Rather, we have eigenvalues collapsing onto an algebraic curve ( Fig. 6), which is inside the unit disk and only tangential to a unit circle at κ = 0 (z = 1). This indicates that the effect of the boundaries is thermodynamically nontrivial. Note that this implies that the thermodynamic limit and the long-time limit do not commute. More specifically, the stationary state is the only asymptotic state for any finite system size n, but with increasing n the time it takes to reach it diverges.

Exact large deviations
The Markovian matrix U describes the ensemble averaged dynamics of a stochastic process. In order to gain further insight into the detailed dynamics of the individual realizations of this stochastic process (also referred to as trajectories) we need to go beyond the ensemble averaged properties. Fortunately there exists an elegant framework [48][49][50] for calculating the fluctuations of time-integrated local (or extensive) observables, ) , (69) in the large time T limit. The ansatz above encodes arbitrary two-site observables that are local (and extensive) both in space and time and are parametrised by functions f x and g x . Note again that the dynamics is defined on a diamond sublattice x + t = 0 (mod 2) (see the discussion around Eq. (1)). This framework is referred to as large deviations (LD), also known as full counting statistics. Before we define it, we just emphasise that the observables of the type (69) depend on the full trajectory realization (s 0 , s 1 , s 2 , . . . , s 2T −1 ), where s t is an n-point configuration on the t-th spacetime saw. For instance, one choice may be the time integrated number of occupied cells It can be shown [48] that for T → ∞ the probability distribution of O T has a large deviation form where ϕ n (x) is called the rate function. This implies that the moment generating function also has the large deviation form The derivatives of the scaled cumulant generating function (SCGF) θ n (σ) at σ = 0 correspond to the cumulants of O T divided by time. The LD functions may be understood intuitively as free-energies of the trajectories and are related to rate functions by a Legendre transform, θ n (σ) = − min x [σx + ϕ n (x)].
In order to compute the SCGF one can show [25,48] that we can deform the propagator U with a tilt or counting field σ, where we introduced diagonal matrices with We then have that θ n (σ) = log λ(σ), where λ(σ) is the eigenvalue of U(σ) with the largest real part. Therefore, by finding this eigenvalue we have access to P T (O) and its cumulants.
The quantity O T explicitly depends on the various values of the trajectory at all times, which makes its exact evaluation a formidable task. Remarkably, we can find an exact solution for its statistics in the long-time limit for arbitrary two-site observables f x (s, s ), g x (s, s ). We again take an ansatz similar to the one we used for the NESS in Subsection 4.1. We write an ansatz for the components of the eigenvector of (72) as, The matrices in the MPA are now site-dependent and we demand that they satisfy the following inhomogeneous bulk algebra, that is generalised to include the tilt operators (74), written for even x: where Z (x) s are site dependent matrices generalising the delimiter matrix S from the homogeneous bulk algebra (42). As before, these matrices admit a 3-dimensional representation, with the precise values of the coefficients v given in Appendix B. We impose the following set of boundary equations for boundary vectors, where the sub-(super-)scripts on U R,L denote the row (column) of the corresponding matrix elements. If these boundary equations are satisfied, the MPAs solve the following eigenvector conditions We then have that λ(σ) = λ R (σ)λ L (σ) is the dominant eigenvalue giving the SCGF. For simplicity we focus on conditional boundary driving (15). The solution to boundary equations (78) consists of boundary vectors (given in [25]) and a consistency condition for the NESS orbital, which takes the form of a quartic polynomial, with and where a n = and Crucially, when σ = 0 the eigenvalue equation (80) reduces to the eigenvalue equation for the NESS orbital (52), as expected. Equation (80) can be solved exactly for arbitrary observables and very large system sizes n. However, we find remarkable universality in the scaling of all observables in the thermodynamic limit n → ∞. As the observables (69) are extensive in the system size n, we have that the limits exist. The SCGF then takes the following simple form where the function ϑ(ς) is such that exp[ϑ(ς)] is the maximum real root of the polynomial Remarkably (85) gives the large system size form of long-time cumulants where · denotes the cumulant. Crucially, for k ≥ 2 we see a divergence for n → ∞ at σ = 0. We can extract explicitly a universal form from (86) of the first few cumulants. From k = 1 we get the average per unit time, 2 (s + s ) and g x = 0) from (69) for n = 10, 30, 500. It approaches the asymptotic form (91) in the thermodynamic limit. The blue curve is the order parameter lim T →∞ O T e −σO T /(T nZ T (σ)) = −θ (σ)/n, which has a firstorder phase transition between phases σ < 0 and σ > 0. The green susceptibility curve θ (σ)/n diverges as ∝ n. The inset is the rate function ϕ(x) with x = O T /(T n). The conditional boundary driving parameters are (α, β, γ, δ) = (1/3, 1/8, 1/2, 2/5).
while from k = 2 we get the corresponding variance per unit time . The divergence of the variance (see Fig. 7) signals a presence of a dynamical phase transition in the statistics of the trajectories. More specifically the scaling function for b + c > 0 has the form, For (b + c) < 0 the asymptotic behaviour is obtained by ς → −ς. We deduce that the SCGF converges to, For details see Ref. [25]. The observables and parameters are the same as in Fig. 7. when b + c > 0. The singularity at σ = 0 corresponds to the first-order phase transition. The behaviour should be contrasted with facilitated models [53] as here there are no fluctuations within each dynamical phase. In order to demonstrate the phase transition we show in Fig. 8, alongside the order parameter as a function of σ, typical active and inactive trajectories (on each side of the transition point σ = 0).

Hydrodynamics of Rule 54
Let us now turn away from the boundary driven setup and consider the dynamics of the system far from the edges, or equivalently, dynamics of an infinite system defined on Z lattice. In this section we start the discussion of the bulk physics by providing the effective hydrodynamic equations that govern the dynamics on the large scales. The starting point is the description of macrostates, from which one extracts the information sufficient to characterize the dynamics on the Euler scale. As an example we provide the parametrisation of the state after the bipartitioning quench and we finish the section by discussing diffusive corrections.

Macroscopic description
We start by introducing n + and n − to denote densities of left and right movers, given as where N ν for ν ∈ {±} are the numbers of particles of both types on the lattice of length n (cf. (33)). To account for the entropic contribution to the macrostate, we introduce n tot ν as the total density of "slots" that can be occupied by particles of type ν ∈ {±}. In the context of thermodynamic Bethe ansatz (TBA) [4,54] (see [55] for a pedagogical review), n tot ν corresponds to the total density (i.e. the sum of the densities of particles and holes). Classically we interpret it as the density of the effective free space in which the particles with a finite width can move [56,57]. As a result of the interactions between particles, the total densities n tot ν exhibit a nontrivial dependence on densities n ± described by the following relation, This identity follows from the TBA description of the quantum generalisation of the model (see [32] for the details). Intuitively, we can motivate it by realizing that two solitons of the same kind are either separated at least by two sites, or there is an oppositely moving soliton in between them, as is demonstrated by the following two diagrams, .
Solitons of the same kind therefore behave as hard rods with a finite length (hence the term −n ν ), but due to the attractive interaction with the other species (i.e. the shift after scattering is −2 sites) their total density increases and we obtain the term +n −ν . This is consistent with the effective description in terms of the flea-gas [57], where the interaction between the particles of the same kind corresponds to the positive jump, while the scattering of oppositely-moving solitons induces the negative jump. We restrict the discussion to the simplest class of generalised Gibbs ensembles (GGE), that are characterized by two chemical potentials, µ + and µ − , corresponding to densities of right and left movers, so that the probability of a configuration s is given by This GGE is equivalent to the state (29) given by the PSA with periodic boundary conditions (see Section 3.3 for the details), and can be alternatively given in terms of the MPS (38) if one replaces the left and right boundary vectors with the trace (see the discussion in Section 7.1). In the large system size limit, the partition function is expressed in terms of the following phase integral involving densities of particles, where the configurational entropy density s[n + , n − ] is given by, In the TBA description of the quantum model, this form of entropy density is precisely the Yang-Yang entropy [32]. From the classical point of view, this form can at first sight be surprising as one could naively expect the second term to vanish [56]. However, the particles in the model do not behave exactly as hard rods, which makes the entropy density more complicated. One can also verify the validity of (97) independently, without relying on the TBA description of the quantum model, by studying the partition sum in finite systems (see Appendix C.2 for details). In the thermodynamic limit the expectation values of densities are given by evaluating this integral using the saddlepoint method, which gives the following relation between n ν and µ ν n tot ν − n ν n ν = e ν , ν = µ ν + log where ν can be interpreted as a single-particle (quasi-)energy [54]. Additionally, for later convenience, we introduce filling functions ϑ ν as the ratios between the density and the total density, Note that the GGE is completely determined either by the pair of chemical potentials {µ + , µ − }, or by the pair of filling functions {ϑ + , ϑ − }.
Excitations on top of these macrostates do not simply move with velocities ±1, but rather they slow down due to the interactions with other particles, and their velocities get renormalized into There are many equivalent ways of finding this result: one can derive it from the definition of the dressed velocity [12,13] by using the TBA description put forward in [32]; one can obtain it by using the appropriate result for hard rods [56,58], or an appropriate soliton gas [57,[59][60][61]; or alternatively the dressed velocity can be identified in terms of expectation values in the GGE [31]. Additionally, the simplicity of the result allows us to obtain it through an elementary argument introduced in [31]: let us imagine that we are tracing a right-moving soliton that starts at time t = 0 at position x = 0 and at time t it arrives to the position x. The soliton moves with the bare velocity 1 and gets displaced for two sites every time it scatters, therefore t − x = 2m s , where m s Figure 9. A trajectory of the slowest right-mover. The red-bordered sites denote the right-mover that undergoes the largest possible number of scatterings, which slow it down to the velocity 1 3 .
is the number of left-movers encountered by the tagged particle. The left-movers are moving with the dressed velocity v − , therefore the tagged soliton will encounter all the left-movers that were at time 0 between 0 and x − v − t, which by identifying v + = x/t gives us, The second relation follows from an analogous procedure involving a tagged left-mover. Note that an additional factor of 1 2 comes from the definition of left/right densities (92). This set of equations provides exactly the dressed velocities given by (100).
For any choice of the parameters {µ + , µ − }, the velocities v ν are restricted to The velocity cannot be larger than 1 (or smaller than −1), since the interactions can only slow the particles down. Alternatively, we can also explain the upper bound by the fact that the dynamics is given in terms of local mutually-commuting operators (cf. (8) and (9)), which restricts the correlation spreading to the causal light-cone spreading with the speed 1 (see the upcoming discussion in Section 7.4). The lower bound is at first sight less clear, but can be easily understood when we recall that two consecutive scatterings of a given soliton should be at least 3 time-steps apart, and therefore its effective speed cannot drop to 0. An example of a configuration corresponding to the soliton moving with the smallest velocity is shown in Fig. 9.

Example: Inhomogeneous quench
Now we have all the necessary ingredients to obtain predictions for the dynamics on the Euler scale. We consider an example of the bipartitioning protocol, where at time t = 0, the left half of the system is prepared in the state determined by (ϑ L + , ϑ L − ), while the state in the right half is (ϑ R + , ϑ R − ). At long times, the state is locally described by a GGE parametrised by ϑ ±,ζ that slowly changes with the ray ζ = x/t, The limiting values far from the junction are given by the left and right GGEs, while for the intermediate values of ζ, the filling functions satisfy the following consistency condition [12,13], Here v ν (ζ) = ν 1+2ϑ −ν,ζ denotes the dressed velocity in the state (ϑ +,ζ , ϑ −,ζ ), as given by Equation (100). The consistency condition describes the profile consisting of two steps, determined by the velocity of left-movers in the state on the right and the velocity of right-movers on the left,
The simplicity of the TBA for RCA54 [32] allows us to find the following explicit form of the coefficients D ν 1 ,ν 2 , Here we obtain an additional factor of 2 w.r.t. the expressions introduced in [31,32,62], due to the convention used in the definition of particle densities (92). The diagonal coefficients D ν,ν can be understood as the variance of the position of the particle of type ν due to the fluctuations in the particle density in the quasi-stationary state, and agree with the prediction of [31,32]. For completeness, we conclude the section by reporting the coefficients A ν 1 ,ν 2 , B ν 1 ,ν 2 , that characterize the nonlinear term, as given by [62], Even thought the expressions are very explicit, we are not able to independently verify the validity of the quadratic term in the Navier-Stokes equation (107), since it is hard to solve it to obtain predictions in this limit.

Space-time duality: time-states and space dynamics
The starting point of this section are multi-point correlation functions of local one-site observables at the same spatial coordinate but at different times, while the system is assumed to be in a stationary state p. Such correlation functions can be interpreted as expectation values of observables that are ultra-local in space, but extended in time, and, intuitively, one can imagine to introduce a corresponding probability vector q (to be defined precisely below), that contains expectation values of all such observables. The construction of an explicit representation of q can be schematically understood as a problem of keeping track of solitons that pass through the origin, as shown in Fig. 10. Generically, there is no guarantee that this is feasible. However, in the case of GGEs introduced in the previous section, the solitons that reach the central site are uncorrelated, which implies an efficient MPA representation of q in terms of 3 × 3 matrices constructed in Subsection 7.2.
The probability distribution q has short-range correlations, therefore a natural question that arises is whether it can be interpreted as a fixed-point (steady state) of another dynamical system, which corresponds to exchanging the roles of space and time in RCA54. In Subsection 7.3 we show that the evolution in space can indeed be given in terms of local and deterministic maps, however, their support is bigger than that of the time-evolution. Equivalently, space evolution can be formulated as a composition of non-deterministic 3-site maps and projectors to a subspace of configurations, as we show in Subsection 7.4, where we also introduce a convenient tensor-network representation of dynamics that provides an algebraic interpretation of the MPA encoding multi-time correlation functions. We finish the section by a short discussion of the physics of one-site correlations.

Stationary states of the closed system
We start by recalling that the NESS introduced in Sections 3 and 4 in the absence of boundary driving (and by assuming periodic boundaries) reduces to a GGE from the previous section -see the discussion in Section 3.3. We will therefore use the parametrisation of this class of stationary states in terms of the matrices W(ξ, ω), W (ξ, ω) introduced in (42). In particular, for a periodic system of even size n, where the partition function Z n is expressed in terms of the transfer matrix T We recall that the chemical potentials associated to the densities of left and right movers are given by µ + = − log ξ and µ − = − log ω.
In this section we are interested in the large system-size limit, therefore we introduce the following graphical notation for the asymptotic (i.e. n → ∞) probability of finite Figure 10.
Schematic picture of the setup. We are interested in probability distribution q of configurations (b 0 , b 1 , . . . , b k ) at the center of the chain x = 0 and different times t, while the system is initialized in the equilibrium state p. Intuitively, this can be understood as keeping track of particles that pass through the origin. Since they are interacting with each other, this is generally a nontrivial task and the computational complexity of q might still be exponentially growing with time. In our case, however, the statistical independence of solitons implies an efficient MPA representation of q. block configurations, where the transfer matrix T = (W 0 +W 1 )(W 0 +W 1 ) is obtained from T by the exchange of parameters ξ, ω. Analogously, the probabilities of configurations with odd block length can be obtained by summing over the value of the last bit, These asymptotic probabilities can be formulated in a more explicit and convenient form by diagonalizing T and T . However, the precise form is not essential for the upcoming discussion, and one can find it in Appendix C.1. These asymptotic probabilities exhibit a nontrivial factorization property: the conditional probability of observing k-th bit to be in the state s k , given the configuration on previous k − 1 sites (s 1 , s 2 , . . . , s k−1 ) only depends on the last three bits. Similarly, the conditional probability of s 1 given the configuration (s 2 , s 3 , . . . , s k ) only depends on (s 1 , s 2 , s 3 ). Explicilty, the following holds, where we remark that this is satisfied independently of parity of the first site and the length k (i.e. it also holds when all the configurations are flipped upside-down). Note that this property, although straightforward to prove (the proof is provided in [27,41]), is nontrivial and for example does not hold for the stationary state of a related cellular automaton [63] (rule 201; cf. Section 9) with otherwise similar properties. A direct physical consequence of the factorization of asymptotic conditional probabilities is statistical independence of solitons: in the stationary state, the probability of finding a left-(or right-) moving soliton is independent of the positions of other solitons, as long as there is no left-mover (or right-mover) on the neighbouring sites. The conditional probability p l of observing a left-mover, if we know that the neighbouring left ray does not contain a left-moving soliton is therefore a well-defined quantity and can be expressed in terms of asymptotic probabilities as follows, where λ is the leading eigenvalue of T = (W 0 + W 1 )(W 0 + W 1 ) (see [27] for the details). Analogously, the conditional probability of finding a right-mover at a given site, if the neighbouring ray on the right is empty, is given by Note that the set of probabilities (p l , p r ) provides an equivalent parametrisation of the steady state, since the mapping (ξ, ω) → (p l , p r ) can be inverted, Comparing these expressions with the TBA equations (98) and (99) we realize that p l and p r coincide with the filling functions ϑ + and ϑ − , respectively.

Time-states
We define time-states as probability distributions of time-configurations, i.e. configurations observed in time (as schematically shown in Fig. 11), under the assumption that the system is in a stationary state p. Explicitly, the components of the time-state q, where s t k is the value at site k of the configuration s = (s −m , s −m+1 , . . . s m−1 ) propagated for t time steps started from initial data s t k = s k , and p s is the asymptotic stationary probability (in our case given by Equations (112) and (113)). Example of time evolution. The red strip denotes a configuration of bits at a fixed time s t , while the orange strip represents an example of a timeconfiguration, i.e. a set of bits observed at the same position in consecutive time-steps. We assume s t to be distributed according to a stationary probability distribution p, while the probabilities of time-configurations s x are given by the corresponding timestate (probability distribution) q.
To understand why for RCA54 the probability distribution q can be explicitly obtained, we imagine the following setup (schematically illustrated in Fig. 10): at time t = 0 we pick a random configuration, distributed according to the stationary probability distribution p, and let it evolve in time, while keeping track only of the bit at position x = 0 or x = 1 (depending on the parity of the time-step, cf. Fig. 11). Generically, the probability of observing a soliton at time t depends on the full observed time-configuration. In our case, however, we recall that the solitons in the underlying stationary state p are statistically independent. This implies that at any time, the probability of a left (or right) mover reaching the observed site, is constant and given by p l (or p r ), as long as there was no soliton in the previous time-step.
Explicitly, the conditional probability of observing a time-configuration (b 0 , b 1 , . . . , b k−1 , b k ) given the previous configuration (b 0 , . . . , b k−1 ) therefore depends only on the last 4 bits, which implies a set of relations analogous to (114), where f, f are two yet unknown functions of the last four bits. Conditional probabilities f ( ) can be straightforwardly determined by observing that short 3-site time-configurations can be divided into three groups.
(i) A configuration can be inaccessible, i.e. does not appear in a time-configuration s x corresponding to any initial configuration s t . To identify them, we only have to note that the following 4 configurations never appear in the local time-evolution rules (3), .
The probability of obtaining such configurations is 0, which means that the conditional probabilities f ( ) (0, 1, 0, s) and f ( ) (1, 1, 1, s) Here the red arrows show the position in which the previously untracked (left or right moving) particle should appear for the new bit in the configuration to be equal to 1.
This exhausts all the possible conditional probabilities, which allows us to express the full probability distribution q as where the last term in the product is f or f for even or odd m, respectively. To completely determine q, we have to take into account the stationarity of the underlying state p, which implies the following consistency condition for any length m, This finally provides (up to normalization) the explicit values of q b 0 b 1 b 2 and we obtain an explicit representation of the time-state in a form similar to the patch-state ansatz introduced in Section 3.2.
The time-state q can be equivalently represented in an MPA form, by encoding the appropriate factors f ( ) in the 3 × 3 matrices A ( ) s , and defining the appropriate boundary vectors L| and |R , Using these definitions, the full time-state can be succinctly summarized as, Note that the right boundary vector |R is the right eigenvector corresponding to the eigenvalue 1 of both A 0 + A 1 and A 0 + A 1 , which immediately gives us access to the normalization of the state,

Space evolution
The discussion in the previous subsection brings us to a related question. Let us imagine we know that for a given trajectory the time-configuration at position x is equal to s x . Is it possible to deterministically map it into the neighbouring configuration s x+1 (schematically represented in Fig. 12)? In other words: we wish to understand whether we can define a local and deterministic model in the space-direction, with respect to which the class of time-states introduced above is stationary. Generally, there is no guarantee that the space evolution can be expressed as a composition of local deterministic maps. In our case, however, we expect this to be possible due to the solitonic description of the model: the dynamics in space can be again understood as solitons moving in one of the two directions with a fixed velocity 1, and undergoing pairwise scattering (see e.g. Fig. 11). The main difference with respect to the usual dynamics in the time direction, is the nature of the interaction, which in this case temporarily speeds the particles up, rather than slowing them down. We therefore expect that it is possible to find a local deterministic map that gives the new value for a bit in the middle of the configuration of 2r + 1 sites, Figure 12. Schematic representation of the mapping between time-configurations. Evolution in space is analogous to time-evolution (see e.g. Fig. 1): at each point x, the bits are updated in rows with the same parity as x, while the others are unchanged.
Note that χ given by (2) is of this form with r = 1. At the moment it is not clear what the support for the local space-evolution maps should be, but by examining the time-evolution rules (3), we quickly realize that the support should be larger than 3 (i.e. r ≥ 2), since e.g. the last two diagrams suggest two different updated values corresponding to the time-configuration (0, 1, 1). Nonetheless, three-site timeconfigurations provide a starting point for construction of space-evolution maps with larger support. Again we realize that three-site configurations can be divided into three types.
(i) As we already noted earlier, time-configurations cannot contain substrings (0, 1, 0) and (1, 1, 1), as these are forbidden by the dynamic rules. The space map φ is therefore defined on the reduced space consisting only of allowed configurations.
(iii) In the case of (0, 1, 1) and (1, 1, 0) one needs more information to be able to find the updated value of the central bit, since in both cases the transitions to 0 and to 1 are possible. However, it suffices to extend the configurations for two sites to top and bottom to find a deterministic update. To demonstrate it, let us first focus on (0, 1, 1). By avoiding the forbidden configurations, there are only two ways to extend the time-configuration for two-sites to the bottom, (0, 1, 1, 0, 0) and (0, 1, 1, 0, 1), .
Now we are able to find a unique way to update these 5-site configurations, using the rules we already have. For either one of these configurations the update rules can be applied to the bottom three sites, after which we try to set the new value at the top to either 0 and 1 and in both cases, only one of them obeys the restriction to allowed configurations, update update , update update . (136) Thus we see that to find the updated value of the middle bit in the subconfiguration (0, 1, 1), it is sufficient to extend the subconfiguration two sites backwards in time.
Similarly, for (1, 1, 0) we have to check the bits to sites forward in time. The rules for all such subconfigurations can be summarized by the following diagrams, where grey squares correspond to the sites that do not influence the updated value (denoted by red-bordered squares).
Together with the restriction to allowed configurations, diagrams from Equations (134) and (137) completely determine the map φ, which has in our case support 7 (i.e. r = 3) and can be summarized as follows, Note that φ(s −3 , s −2 , s −1 , s 0 , s 1 , s 2 , s 3 ) does not depend explicitly on the values s −2 , s 2 , which means that all the maps applied in the same step (see Eq. (133) and Fig. 12) commute. The dynamics in space can be therefore understood as a composition of strictly local maps, that are in each step applied to all 7-site subconfigurations centred around the sites labeled by temporal indices of the same parity, and the order in which they are applied is not important. From this point of view, space-dynamics is, apart from the larger support of local maps, defined in perfect analogy to the usual evolution in the time-direction.

Circuit representation of dynamics
The staggered structure of time evolution, where each time-step is given in terms of mutually commuting local updates, implies a natural representation in terms of a (classical) tensor-network. The motivation for this comes from recent advances in constructing (typically nonintegrable) solvable models of many-body nonequilibrium dynamics in terms of quantum circuits, with notable examples being random [64][65][66][67] and dual-unitary circuits [68][69][70][71][72][73]. These models constitute a family of analytically tractable quantum many-body systems, for which many previously not-accessible outof-equilibrium quantities can be obtained exactly. As we show below, RCA54 admits a convenient representation in terms of local gates, which allows us to find an equivalent diagrammatic procedure to recover the results from the previous subsections. We start by defining two tensors represented by big and small circles, where each leg can be either in a state s = 0 or s = 1, The big circle encodes the deterministic update (2), while the small circle forces each one of the legs to be in the same state. Using this graphical convention, the local time-evolution operator (10) can be represented as The definition of the smaller tensors immediately implies that two circles sharing a leg can be combined into one with more legs. In particular, we note the following identity, This allows us to combine all the local operators applied at the same time-step into one horizontal line of small and big balls positioned at sites with the opposite parity, U e = · · · · · · = · · · · · · , (142) and U o takes a similar form with exchanged roles of the two tensors. Interchangeably applying U e and U o we obtain a network consisting of small and big tensors, positioned at sites with the sum of the two coordinates x + t ≡ 0 (mod 2) and x + t ≡ 1 (mod 2), respectively, as is diagrammatically shown in Fig. 13. tensorÛ , defined as, OperatorsÛ e andÛ o , describing the full one-step evolution in space, are in analogy to U e and U o given in terms of products ofÛ centred around the sites on the dual lattice with even and odd parity respectively, as is diagrammatically shown in Fig. 13. One quickly realizes thatÛ is not deterministic, which is expected, since we know that the deterministic map has the support 7. To understand how this happens, we define the following 3-site projector to the space of allowed configurations, We note that the local three-site dual operator projects to the same subspace as P , and additionally, P andÛ commute if they overlap for at most one site, while they do not commute if the overlap is 2, This allows us to show that the evolution in space can be equivalently given in terms of the operators reduced to the space of allowed configurations, where P e/o are defined as products of P on even/odd sites (in analogy toÛ e/o ). To get from the first to the second diagram, we use the relation (145), and to obtain the third diagram we note that all projectors P commute amongst themselves. The operatorsŨ e ,Ũ o are precisely the space-evolution maps on the reduced subspace, and from the previous subsection we know that they can be expressed as a composition of 7-site deterministic maps. One can show this independently by expressing the 7site deterministic gates in terms of the non-deterministic gateÛ squeezed between the relevant projectors, and then prove that the two descriptions are equivalent by utilizing appropriate local algebraic relations fulfilled by the gates and projectors. The full construction with the proof is reported in Ref. [28].

Tensor-network representation of multi-time correlation functions
To demonstrate the usefulness of the circuit formulation of RCA54, we return to the multi-time correlation functions at the same point and show an independent derivation of timestates. For simplicity we will only consider the maximum-entropy state, but a similar analysis can be done for a stationary state of the form introduced in Subsection 7.1. Before starting, we recall that by definition the expectation value of an observable a in a probability distribution on n sites p is a sum over all the configurations of the products of the probability and the value of the observable in the configurations, where we introduced the unnormalized one-site maximum-entropy state ω and A is a 2 n ×2 n diagonal matrix with the values a(s). Note that the maximum entropy state on n sites corresponds to p = 2 −n ω ⊗n . This allows us to introduce the multi-time correlation function of one-site observables a 1 , a 2 , . . . , a k at times t = 0, 1, . . . k − 1 as, where represents a one-site observable and = ω. For simplicity we will assume that k ≡ 1 (mod 2), but the same can be repeated for even k. Note that these correlation functions can be directly related to time-states q as To reduce the diagrammatic representation of C a 1 ,...,a k , we first note that the local time-evolution operator U is deterministic and as a consequence the maximum-entropy state is invariant under it (as well as under its transpose), Therefore the circuit (149) can be simplified by removing gates connected to ω, which reduces it to a tilted square shape. Additionally, since the observables are diagonal, they commute with the small tensor, and we can recast the circuit in terms of dual gatesÛ as So far we have not taken into account any special property of RCA54 and the computational complexity of such an object would generally still grow exponentially. However, we know that in our case the complexity is constant, therefore there must be some simplification occurring.
If the gatesÛ were deterministic, a relation analogous to (151) would be satisfied in the space direction and the correlation functions would trivialize. This does not hold precisely in our case, but the fact that the dual evolution is deterministic on the reduced subspace suggests that the gatesÛ have some additional structure. In particular, we find that the following two few-site algebraic relations are fulfilled, and sinceÛ T =Û , an analogous set of left-right flipped identities holds. These relations are used, together with the projector identity (145), to remove dual gates from the diagram (153) layer by layer, starting at the edges and working our way inwards. We are left with only two layers of gates, one on each side of the column of observables, at which point the algebraic relations (154) can no longer be applied. Note that here we implicitly used the fact that observables are diagonal and therefore commute with projectors P . By definition (cf. (150)) the fully simplified tensor network (155) gives the following diagrammatic representation of components of the time-state q corresponding to the maximum entropy state, The diagram immediately implies that the Schmidt rank does not grow with time and can be bounded from above by 4. The equivalence between this object and the corresponding time-state (130) can be straightforwardly shown by explicitly extracting 4 × 4 matrices and boundary vectors from (156). The details of the proof are reported in Ref. [28], together with the generalisation of the circuit representation of the multi-time correlations to the full class of Gibbs-like stationary states p. Furthermore, analogous treatment can be straightforwardly recast for the quantum version of the model [35][36][37]. The main conceptual difference between the quantum and classical case is the size of the vector space, which is in the quantum case doubled (squared in dimensionality), as we are considering density matrices rather than probability distributions. Despite that, a similar set of few-site algebraic relations can be formulated, which gives access to objects analogous to (156).

Correlation functions at one site
The MPA representation of q gives us information about all multi-point correlation functions. In particular, one can obtain the full probability distribution of gaps between the consecutive observed particles, which in related stochastic models shows nontrivial phenomenology [53,74,75]. In our case, however, this behaviour is not reproduced, since the solitons passing through the origin are uncorrelated (see [27] for details).
Another quantity of interest is the one-site dynamic density-density correlation function where ρ denotes the local density of full sites at position x = 0. The correlation function can be easily expressed in terms of the MPS (130) as, where we for simplicity assume even t. Since the matrices A ( ) s are finite-dimensional, the above matrix product can be easily evaluated and shown to be equal to where Λ 1,2 are subleading eigenvalues of the transfer matrix (A 0 + A 1 )(A 0 + A 1 ), The correlation function decays exponentially for all values of parameters, and furthermore, one can show [35,36] that the correlations exponentially decay also when we move away from x = 0 as long as |x/t| < 1 3 . This is compatible with the hydrodynamic picture, according to which the correlations move with velocity larger than 1 3 (cf. (102)), and therefore we cannot observe power-law decay on rays with |x/t| < 1 3 .

Time-dependent matrix product ansatz
In the previous section we were discussing evolution of an observable restricted to one spatial coordinate. Here we go one step further and discuss the full time-evolution of local one-site observables, which is the classical analogue of the Heisenberg picture timeevolution. In generic systems, this is a hard problem and the computational complexity typically increases exponentially with time. However, in the case of Rule 54 the problem reduces significantly and the complexity of time-evolution grows quadratically with time. This is the consequence of the fact that the deterministic time-evolution reduces to the problem of back-tracking positions of solitons, which can be due to the simplicity of the interactions done efficiently. Most of the section is devoted to the construction of the time-dependent matrixproduct state that encodes this procedure, which allows us to find exact density profile after an inhomogeneous quench and the density-density spatio-temporal correlation function in the maximum entropy state. We conclude by discussing the quantum generalisation to non-diagonal observables, which provides a strict bound on operatorspreading in the quantum version of the model.

Dynamics of classical observables
We assume the finite chain with periodic boundaries of even length n, and the sites are labelled by integers in the interval [−n/2 + 1, n/2], so that the site 0 is at the middle of the chain. The length of half of the chain n/2 is assumed to be strictly larger than any finite time t we are considering. Here we use U(t) to denote time-evolution for the first t time-steps, i.e.
and in the first time-step, even sites get updated.
As we already noted in the previous section, classical observables can be understood as diagonal operators acting on the 2 n -dimensional space of (macroscopic) states, so that the classical analogue of the Heisenberg picture time-evolution is given by which in our case reduces to a(t + 1) = U e a(t)U e , t ≡ 0 (mod 2), In the following we will consider the simplest local observable: local density in the centre of the chain ρ = ρ 0 , where ρ x denotes the density at the position x, Here we introduced the shorthand symbol |s s| x , s = 0, 1, for the projector to one of the local physical states at site x. Additionally, we assume the following notation for the projector to the configuration Note that ρ x (t) can be obtained from ρ(t) by a lattice spatial shift operator η, Furthermore, ρ x (t) also gives immediate access to the time-evolved local density of empty sitesρ x = |0 0| x , through the relationρ x (t) = 1 − ρ x (t). Therefore, assuming that we know ρ(t), we are in principle able to express time evolution of any local diagonal operator, since time evolution of observables is a homomorphism, (AB)(t) = A(t)B(t). Time evolution is deterministic and local (see e.g. (151)), which implies that at time t the time-evolved observable acts nontrivially on the section of the chain between −t and t, where c s (t) are the appropriate coefficients in the basis expansion. By the definition of time evolution (10) each one of these observables will be mapped into observables with a (two sites) larger support, where the choice between U e and U o depends on the parity of the time step (cf. (163)). This immediately implies that there are exactly 4 t different configurations s, referred to as accessible configurations, of length 2t + 1, for which c s (t) = 1, while the other 4 t coefficients are zero. Furthermore, accessible configurations are precisely the ones for which the site at the origin was full at time t = 0, which means that the configuration s contains a soliton that was at time t = 0 at the site x = 0. The problem of time-evolution of local density can be therefore mapped onto the equivalent problem of identifying solitons in the configuration and backtracking their position in time. By itself, this realization is not very deep and in principle the complexity of timeevolution could still grow exponentially with time t. However, in our case, the dynamics of solitons is very simple, which significantly reduces the complexity of the problem. In particular, for any soliton found in the configuration s = (s −t , s −t+1 , . . . , s t ) it is possible to determine whether or not at time t = 0 it was at the origin by counting the number of scatterings it was involved in. An example of an accessible configuration is shown in Fig. 14. It contains several solitons, amongst which there is a left-mover (denoted by the red colour), that was at time t = 0 at the position x = 0. To show that the red soliton passed through the origin, one has to identify three right-moving solitons it scattered with, which are all contained in the part of the configuration to the right of it. However, one should be careful not to count too many right movers, as e.g. the orange soliton never scattered with the red one, which is achieved by keeping track of the left-movers positioned to the right of the red soliton.
To make this simple idea more precise, we start by splitting the coefficient c s (t) into the sum of two terms c L s (t) and c R s (t), which give 1 when the soliton from the origin is moving in the left and the right direction respectively. As we will argue later, both contributions can be efficiently expressed in terms of products of matrices where X L/R s , Y L/R s ∈ End(V), s ∈ {0, 1} are linear operators over the infinite dimensional auxiliary Hilbert space V, V = span |c, w, n, a | c, w ∈ N 0 , n ∈ {0, 1, 2}, a ∈ {0, 1} , and L L (t) , R L , R R (t) , R R boundary vectors from the auxiliary space. Since the boundary vectors change with time t, we refer to the expression (169) as the timedependent matrix-product ansatz (tMPA). To compactly express the matrices, we define the ladder operators c +/− , w +/− that change the value of c and w by one, c + = c,w,n,a |c + 1, w, n, a c, w, n, a| , c − = c + T , w + = c,w,n,a |c, w + 1, n, a c, w, n, a| , and the projectors c c 1 c 2 = w,n,a |c 1 , w, n, a c 2 , w, n, a| , w w 1 w 2 = c,n,a |c, w 1 , n, a c, w 2 , n, a| , n n 1 n 2 = c,w,a |c, w, n 1 , a c, w, n 2 , a| , a a 1 a 2 = c,w,n |c, w, n, a 1 c, w, n, a 2 | .
We will elaborate on the physical interpretation of the auxiliary space in the next subsection. The operators X L s , Y L s can be expressed as 3 × 3 matrices acting on the space spanned by {|n } 2 n=0 as a 00 c − w + + a 11 w + 0 0 a 00 c − w + + a 01 w + + a 11 c + w + 2 0 0 a 00 c − w + + a 11 w + 0 0 while the time-dependent boundary vectors are given by The contribution c R s (t) can be obtained from c L s (t) by taking a transpose, and adding additional boundary terms (see the discussion at the end of this section for details), The structure of these 3 × 3 matrices is very similar to the representation of the cubic algebra introduced in Section 4.1. Indeed, operators (42)), where the entries are no longer scalars, but rather more general operators acting on an infinite-dimensional vector space.
Before discussing the construction of tMPAs for c L s (t) and c R s (t), let us remark that even though the operators X L/R s , Y L/R s are acting on an infinite-dimensional Hilbert space, we can for every finite t replace them by finite matrices, with the dimension that scales as O(t 2 ). In particular, the possible values of c and w in definitions of ladder operators (171) and projectors (172) can be restricted to 0 ≤ c, w ≤ t + 1. Note that this is an exceptional property, since typically one expects the complexity to grow exponentially with time.

The contribution of left-movers
We start the construction of the matrix-product representation (169) by discussing the first term c L s (t). Afterwards we will show how to adapt the result to obtain also the contribution from the right-movers.
Let us consider a configuration (s −t , s −t+1 , . . . , s t ) = s ∈ Z 2t+1 2 . To figure out whether or not c L s (t) = 1, we imagine starting at the left-most site of the configuration and moving towards the right, reading the configuration s site by site. Whenever we encounter a left-mover, we designate it the test soliton and, aiming to determine whether or not it originated from the centre, we start counting the solitons on its right. To encode the soliton counting procedure, we introduce four auxiliary degrees of freedom, |c, w, n, a .
(i) The activation bit, a ∈ {0, 1}, tells us whether we are on the left (a = 0) or the right (a = 1) side of the test soliton. If the activation bit is turned off, the state splits into two parts whenever we encounter a left mover. The first part corresponds to a = 0, describing the situation in which the left mover is not the test soliton, while the second part represents the opposite case, with a = 1. In any other situation the activation bit remains unchanged.
(ii) The collision counter, c ≥ 0, represents the number of scatterings the test soliton had to undergo if it passed through the origin. As long as a = 0, the collision counter increases by 1 every two sites. If a = 1, the collision counter decreases by 1 whenever a right-mover that scattered with the test soliton is encountered. When we reach the right edge of the configuration, the value c = 0 tells us that the test soliton passed through the origin, and c = 0 implies the opposite.
(iii) The scattering width, w ≥ 0, keeps track of the number of scatterings of the rightmovers encountered after the test soliton. At the left edge the width w is equal to the time-step t, and every two sites it decreases by 1. Additionally, the width changes as w → w − 1 whenever a left mover on the right side of the test soliton is encountered. All the right-moving solitons that we meet after w drops to 0 could not have scattered with the test soliton.
(iv) The occupation counter, n ∈ {0, 1, 2}, provides additional information about the particle content needed to appropriately change w and c. Explicitly, n = 0 if the current site is empty, n = 1 if the site is full and the left neighbour is empty, and n = 2 if the site and the left neighbour are both occupied.
At the left edge, the collision counter should be 0, and the scattering width is set to be equal to the time-step t. Furthermore, we start with the activation bit equal to 0 (since we have not yet met any test soliton), and the particle content outside of the lightcone does not influence the soliton counting, therefore we can without loss of generality assume n = 0. This implies the following form of the left boundary vector, The right boundary vector should have nonzero overlap with vectors that correspond to a test soliton that passed through the origin, which is given by c = w = 0 and a = 1, while n can be arbitrary, |0, 0, n, 1 .
To construct the matrices that correspond to the procedure summarized above, we consider 3 regimes.
Left side of the test soliton. Before the test soliton is encountered, the scattering width decreases by 1 every two sites and at the same time the collision counter should increase, while n should be appropriately adjusted to keep track of the consecutive full sites. The left action of the matrices on the sector a = 0 is therefore, c, w, n, 0| X L s a 00 = c, w, s · min{n + 1, 2}, 0| , c, w, n, 0| Y L s a 00 = c + 1, w − 1, s · min{n + 1, 2}, 0| .

(178)
Encountering the test soliton. Whenever a left-mover is encountered while a = 0, an additional vector with a = 1 is created. There are 4 configurations corresponding to this situation. The first two are simpler and represent a situation where a left-mover is detected while moving uninterrupted, , which can be summarized with the following matrix elements c, w, 1, 0| X L 1 a 11 = c, w, 2, 0| X L 1 a 11 = c, w, 2, 1| .
The other two configurations describe the soliton encountered while it is scattering, , where the grey squares denote the path of the soliton in the previous time-steps. The corresponding matrix elements are where in the first case we note that one scattering already occurred (hence c − 1), while in the second case the decrease of w is just due to moving to the right. This exhausts all the subconfigurations in which the bit a can flip (from the left) from 0 to 1.
Right side of the test soliton. Let us first assume w > 0. After the test soliton has been chosen, the changes of c and w are summarized as: Explicitly, there are two possible configurations of an uninterrupted right-mover appearing, , which are described by, while the configurations in (179) imply the following, The configurations in (181) describe the two scattering solitons, therefore the corresponding matrix elements contain the decrease of both c and w, All the remaining configurations do not contain any newly detected solitons, , therefore c remains constant and w decreases only due to moving two sites to the right, The right-moving solitons encountered after the width drops to 0 could not scatter with the test soliton, because they were too far. Therefore c should stop decreasing, Combining the matrix elements from Equations (178-189), we finally recover the matrices given in (173). This completes the construction of the tMPA corresponding to c L s (t).

The contribution of right-movers
The tMPA for c R s (t) can be derived in a similar fashion, by simply reversing the directions of all solitons, which is achieved by exchanging the roles of the boundary vectors and by transposing the tMPA matrices. Additionally, we have to exclude all the configurations with a central right-mover that were already captured by c L s (t). These are exactly the configurations where at time t = 0 the soliton in the origin is in the process of scattering, i.e. their trajectory starts with one of the following two diagrams, . ( In the previous case, these would correspond to the configurations with another scattering occurring exactly when w = 0. Therefore to exclude them, we have to increase w in the right boundary vector by 1, while at the same time modify the left boundary vector to allow for w not dropping to 0, R R (t) = L L (t + 1) , L R = R L + 0, 1, 0, 1| + 0, 1, 2, 1| . (191) Additionally, it does not suffice to map X L s T → X R s , Y L s T → Y R s , but we have to further modify them so that the following matrix elements are 0, which gives exactly the form summarized in (175).

Physical applications
The tMPA provides means to study transport in the model. In this subsection, we provide two examples of physically relevant quantities that can be exactly obtained with our solution: density profile after a bipartite quench and a dynamical densitydensity correlation function. We finish the subsection by suggesting a formal extension of the tMPA to richer stationary states.
8.2.1. Inhomogeneous quench Let us consider a particular bipartitioning protocol, where at time t = 0, a half-infinite chain prepared in the maximum-entropy state is joined together with a half-infinite empty chain, so that the initial state is given by, We wish to obtain the time-dependent profile of the particle densityρ(x, t), which is given by the expectation value of ρ x in the time-evolved inhomogeneous state, Note that the negative sign in ρ x (−t) together with the staggering forces a different convention with respect to (166). The locality of time-evolution immediately implies that outside of the light-cone the density profile matches the appropriate boundary values, i.e.ρ(x < −t, t) = 1 2 and ρ(x > t, t) = 0. For intermediate values of x, the expectation value is by definition (see Eq. (148) and the discussion around it) where ω = [1 1] T denotes a (unnormalised) one-site maximum-entropy state. The expectation value is (assuming odd x + t, cf. (194)) conveniently expressed in terms of the tMPA as The matrices X L/R s , Y L/R s are infinitely dimensional, therefore it is not immediately clear whether it is possible to evaluate this expression. However, their sparse structure enables us to find exact form of certain simple matrix elements (see [26] for the details). In particular, one can prove that r(m, t) takes the following form, where θ x denotes the discrete Heaviside function (i.e. θ x≥0 = 1 and θ x<0 = 0).
If m < (2t + 1)/3, only the last term remains and r(m, t) greatly simplifies to give the following density profile for x > −(t + 1)/3, This profile consists of a front moving ballistically to the right with velocity 1 and exponentially suppressed oscillations behind it, as can also be seen in Fig. 15. Additionally, it exhibits an intuitive physical picture: in this section of the light-cone, all the solitons are moving to the right with the dressed velocity 1, since the right half of the lattice was empty at the beginning and there are no left-movers to slow them down. To determine the left-edge of this section, we note that the right-most left-moving solitons are moving with the velocity −1/3 (cf. Fig. 9 and the discussion in Section 6), therefore the full section of the lattice between x = −t/3 and x = t is populated only by right-movers [26]. As we increase m (or equivalently, as we move towards the left edge of the lightcone), we are no longer able to make further simplifications to the expression (197). However, one can easily find an asymptotic expression that describesρ(x, t) on large scale. Using the Stirling formula it is possible to show that the density profile is well approximated by a shifted error function centred around x/t = −1/2, with the width that scales as √ t, On the rays with the constant ζ = x/t, the density profile matches the Euler scale prediction (106): the profile consists of two steps, with their velocities equal to − 1 2 and 1. In the central part, the hydrodynamic prediction is given by the expectation value of density in the state determined by (ϑ + , ϑ − ) = ( 1 2 , 0), and can be shown to be equal to 1 3 using the MPA representation of GGEs (see Appendix C). Furthermore, the diffusive broadening around the left step is compatible with D −− = 1 16 , which is consistent with the expression (108) in the case ϑ + = 1 2 . The right-step does not exhibit any diffusive scaling, since the term D ++ vanishes for ϑ − = 0. The oscillations around 1 3 decay exponentially with the distance from the step and are therefore not detectable on the hydrodynamic scale.

Dynamical structure factor
We proceed to the second example: the dynamic density-density correlation function, C(x, t), evaluated in the maximum-entropy state, where we introducedC(x, t) to denote the rescaled expectation value of the product of densities. Analogously to r(m, t) defined before,C(x, t) can be expressed as a simple sum of matrix products corresponding to the left and right tMPAs, where we assumed x + t ≡ 0 (mod 2). In the opposite case, we take the advantage of the identity C(x, t)| x+t≡1 (mod 2) = C(x, t − 1), and thus reduce it to the previous case. It is again possible to evaluate these sums, and we obtain the following expression [26], where for simplicity of notation we assume n<0 k = 0 for any k, n ∈ Z. The correlation profile can be again split into two qualitatively different parts. If |x| < (t + 1)/3, the spatial dependence disappears (apart from the staggering) and the correlations reduce to In this regime, correlations decay exponentially with time, |C(0, t)| ∼ 2 −t/2 , and the above expression exactly matches the infinite temperature result obtained from multi-time correlation functions in Section 7. Furthermore, in this regime C(x, t) is homogeneous in space and exponentially suppressed also when one considers the underlying stationary states to belong to the richer class of two-species GGEs considered in the previous sections, as a consequence of the minimal velocity of excitations being equal to 1/3 (cf. (102)). This can be proven by generalising the circuit approach outlined in Section 7 [35,36]. If |x|/t > 1/3, we are left with a binomial sum that cannot be further simplified, therefore it again makes sense to find the asymptotic shape of the correlation profile. As is shown in Fig. 16, the profile consists of two ballistically moving peaks with velocities ± 1 2 that spread diffusively. Indeed, one can straightforwardly show The result nicely matches the hydrodynamics: the peaks move with dressed velocities v ± = ± 1 2 and their diffusive broadening is governed by D νν = 1 16 , ν ∈ {+, −}. Both quantities are compatible with their predictions given by (100) and (108) in the maximum-entropy state. s , allows us to simply expand the tMPA to include the right probabilities by a simple element-wise multiplication of both sets of matrices, and thus generalise the approach to the class of GGEs given by (95). To make this point more precise, let us consider the dynamical correlation function in a state p given by the pair of parameters (ξ, ω) (corresponding to µ + = − log ξ, µ − = − log ω in (95)), which is by definition expressed as where the state-dependent coefficient c * s (t) is defined as the product of the tMPA coefficient c s (t) and the probability of finite configuration (113). Using the MPA formulation of these probabilities outlined in Appendix C.1, the state-dependent coefficient can be rewritten as Since both the matrices describing stationary probabilities and the matrices constituting the tMPA have very similar 3 × 3 block structure (cf. (173), (175), and (42)), we can in analogy to (169) express c * s (t) as a sum of the left and right time-dependent matrix products, by defining the new operators X and similarly, Here denotes the element-wise multiplication, i.e. a 1 a 2 a 3 The generalised expressions (207,208) can be understood simply as attaching a weight ξ or ω to each left or right-mover that we observe in the soliton-counting procedure. When ξ = ω = 1 we recover the original tMPA matrices and vectors (173,174,175). Note that in the case of right-movers we have to use the transpose of the GGE matrices, and additionally the role of parameters is swapped. The explicit form of the new set of operators can be found in [41].
This construction in principle allows us to obtain C(x, t) for the full class of twoparameter GGEs (space and time translation invariant states) and, by appropriately choosing parameters, one can also obtain the tMPA formulation of the time-dependent density profile after a bipartite quench from an arbitrary pair of these GGE states. However, so far no attempt has been done to evaluate these expressions and obtain explicit results. Whether or not this is feasible is still an open question.

Quantum interpretations and operator spreading
Even though we have been discussing RCA54 in the context of classical dynamics, there is nothing preventing us from treating it in the quantum setting. Indeed, since the one-time-step evolution operators U e and U o describe deterministic dynamics, they are also a special case of unitaries, (U o U e ) (U o U e ) † = 1. In this respect the model represents one of the simplest interacting models, in which different measures of operator spreading can be easily studied [30,31,33,34]. In particular, one such quantity is operator space entanglement entropy (OSEE), which measures the complexity of simulability of quantum dynamics in the Heisenberg picture [76][77][78][79]. In the case of RCA54 the bound of the rate of OSEE growth with time can be obtained exactly, by generalising the tMPA introduced above to the time-evolution of quantum (rather than classical) observables [33].
Time-evolution of a quantum one-site observable a can be split into 4 contributions, each one corresponding to a basis operator |s 1 s 2 |, a = a 00 a 01 a 10 a 11 , a(t) = U(t) † aU(t) = a 00 U(t) † |0 0| U(t) + a 01 U(t) † |0 1| U(t) We immediately note that the diagonal elements, , except for the initial site, which is 1 for s and 0 for b. The sites denoted by the red, yellow and orange colour in both cases denote the same three solitons, while in the top part the grey sites denote solitons, whose positions are the same as in the bottom part. To understand how s maps into b, we note that in this example, changing the initial state from 1 to 0 causes an additional right-mover to appear (between the yellow and orange soliton), which causes the displacement of all the solitons that scatter with it.
are given by the tMPA introduced in the previous subsections, while one of the offdiagonal terms can be expressed as the transpose of the other, Therefore, the only remaining thing to do is to find an efficient matrix-product description of σ(t) = U(t) † |0 1| U(t).
Analogously to the diagonal case, at time t the time-evolved operator σ(t) nontrivially acts on the sublattice of length 2t + 1 and we can introduce coefficients By definition, each one of these terms will be at time t + 1 mapped into 4 different operators with support 2(t + 1) + 1, where we use the shorthand notation s x = χ(s x−1 , s x , s x+1 ) and b x = χ(b x−1 , b x , b x+1 ). The deterministic nature of time evolution has direct implications that can be summarised in the next three points.
(i) The coefficient d s,b (t) can be either 0 or 1.  An example illustrating the point (iii) is shown in Fig. 17. In this case, changing the initial bit from 1 to 0 produces an additional right-mover, therefore s can be mapped into b by adding the additional soliton and displacing all the quasi-particles that were affected by it. More generally, it can be shown that changing 1 to 0 in the initial configuration can be always understood as adding or removing 1 or 2 solitons, and in all these cases, there exists a simple deterministic map M t from s to b that accounts for the relevant displacements of quasi-particles. Furthermore, as was demonstrated in Ref. [33], the map M t can be expressed in terms of products of finite-dimensional matrices, with the size that does not change with time t. This immediately implies that the growth of the complexity of time-evolution is given by the size of operators constituting c s (t), which is O(t 2 ). The growth of OSEE is therefore bounded from above by 2 log(t) (up to a constant factor). Note that this bound need not be saturated, and indeed, there are recent indications of OSEE increasing as 1 2 log t [34].

Related exactly solvable models
While, as argued in this review, the RCA54 model allows for a remarkable set of explicit results on nonequilibrium dynamics and statistical mechanics, an important question is if similar can be achieved for other, perhaps closely related models. Similar results have been so far demonstrated in two other (RCA type) deterministic lattice models. Specifically, in Refs. [63,80] the rule 201 RCA has been studied, which can be understood simply as a negation of rule 54, namely with the local rule modified as χ 201 (s, s , s ) = 1 − χ 54 (s, s , s ).
Both the RCA201 and RCA54 can be interpreted as deterministic variants of kinetically constrained (or facilitated ) models, typically studied in the context of glassy dynamics [20,[81][82][83][84]: the middle bit (cell s ) flips whenever the neighbouring bits (cells s, s ) satisfy a particular condition. In the case of RCA54 the flip occurs when at least one of s and s is in the state 1, while in RCA201 the flip happens when both s and s are in the state 0, and as such it provides a deterministic classical analogy of the PXP model [85][86][87][88]. It has been shown [63] that the invariant (equilibrum) states and NESS of the boundary driven setup can again be written in terms of a patch state and matrix product ansatz with very similar constituent matrix algebra as for the rule 54. However, at the moment it is still unclear precisely which other features of RCA54 dynamics remain exactly solvable for the RCA201, which is left as an open question. Similar explicit results are expected to be achievable for the free RCA rule 150, χ 150 (s, s , s ) = s + s + s (mod 2) (or equivalently, the flip occurs when exactly one of s, s is 0 and one 1), which represents a caricature of non-interacting dynamics in 1 + 1 dimensions. Another class of related deterministic discrete nonequilibrium dynamics that has been studied in some detail is the two-species cellular automation describing a synchronous lattice gas of charged point particles with hard-core interaction [89][90][91]. This model fundamentally differs from the RCAs of Ref. [18] (such as rules 54,150,201) in the way the global time-evolution is defined. In the first time-step, the even pairs of neighbouring cells are updated using a two-site update rule, (s t+1 x , s t+1 x+1 ) = φ(s t x , s t x+1 ), s t x ∈ {+, −, 0}, while in the second time-step, the update rule is applied to odd pairs of sites, and the two steps are then periodically repeated. This is analogous to trotterisations of nearest-neighbour interacting models, or to brickwork circuits. The two-species hardcore interacting model specified by a self-invertible bijective mapping, φ(±, ∓) = (±, ∓), and φ(s, s ) = (s , s) otherwise, has been shown to display a coexistence of ballistic and diffusive transport (similar to Rule 54) with explicitly calculable diffusion constant and Drude weight [89,90], explicitly solvable boundary driven NESS exhibiting a nonequilibrium phase transition [91], but with a slightly lower complexity of time evolution compared to RCA54: the Schmidt rank of time-dependent MPA grows as ∝ t [91] rather than t 2 .

Conclusions and perspectives
This review paper provides a coherent overview of a series of explicit results, appearing over the last five years -most of which were co-contributed by the authors, on dynamics and nonequilibrium statistical mechanics of a particular interacting reversible cellular automaton (RCA54). We believe that establishing such a set of explicit results on the connection between microscopic reversible laws of motion and effectively irreversible dynamics of macroscopic states is indispensable for the field whose holy grail is to establish such a connection on a general level [17].
The phenomenology that the model helps to rigorously understand is the emergence of Fick's law of diffusive transport, as a leading-order correction to the ballistic motion. One can thus consider the RCA54 as an exactly solvable example of nonequilibrium universality class of diffusive transport. However, it has been recently discovered, that even without any intrinsic sources of noise (or dissipation), one-dimensional systems can display other types of super-diffusive (but sub-ballistic) transport, for example in the presence of integrability and non-abelian global symmetries, spin (charge) transport seem to follow the scaling of KPZ universality class with dynamical exponent 3/2 (rather than 2 for the diffusive class) [92][93][94][95][96]. It would be an interesting challenge to look for minimal solvable models of such anomalous transport within a class of reversible cellular automata. For instance, one can generalise RCA to multi-colour internal states of nonempty cells such that it reduces exactly to RCA54 if the colours are disregarded. The conserved charges of such automata (some of which appear to be integrable) seem to display superdiffusive transport [97].
A different potentially interesting direction is to study quantum and/or stochastic deformations of Rule 54.
Specifically, one can allow for coherent (quantum superposition) or stochastic mixtures of several local processes which may preserve exact solvability. For instance, replacing the local (3-site) propagator (10) by a nonpermutation matrix which maps (000) → u 00 (000) + u 01 (010), (010) → u 10 (000) + u 11 (010) (while all other processes remain as defined by (10)), where the 2 × 2 mixing matrix u ij is either unitary or stochastic, one finds strong circumstantial evidence of integrability in all such cases [98]. Another kind of quantum deformation of RCA54, introducing a non-trivial quasiparticle dispersion relation, was shown to allow for standard coordinate Bethe ansatz treatment [32].
Perhaps the most urgent question regarding RCA54 is how to precisely fit it into the family of Yang-Baxter solvable models. At the moment it seems that the model can be interpreted as a (in some sense) degenerate limit of a more complicated integrable system. One example of the mapping in this direction is the quantum deformation [32] mentioned above, in the context of which RCA54 is a zero-dispersion limit of a more general Bethe-ansatz solvable model. Another possibility [99] is to identify it as a singular limit of a vertex model [100], by first mapping it to a randomtiling model [101,102]. However, at the moment it is still not clear what is the best way to fit everything together. Nonetheless, interpreting RCA54 as a singular limit of a more generic model, puts it into the same family as other special (but not free) points of interacting models, such as the q → ∞ limit of the q-boson model studied in [103,104], or the model considered in [105][106][107], with which it shares some common features. It would be interesting to understand this connection in more detail.
Another interesting perspective for future research is studying non-stationary dynamics in the RCA54. More specifically, this would involve the existence of extensive quantities A (sums of local terms) satisfying a Floquet dynamical symmetry [42,[108][109][110] condition UA = e iω A, where ω = 0 is real, which would imply persistent oscillations and absence of relaxation.

Appendix B. Explicit form of the representation for the large deviation MPA
The coefficients in (77) are given by the following explicit two-parameter (ρ, κ) form (ρ and κ play the role of spectral parameters, see Ref. [25] for details), and recast the requirement (C.9) in terms of identities that have to be satisfied for coefficients α n/2 (x, y). It can be then verified that the solution takes the form α k (x, y) = k(k + x + y) (k + x − y)(k − x + y) k + x − y y k − x + y x , (C.11) where we implicitly assume α k (x, y) = 0 for k ≤ ±(x − y). Asymptotically, the combinatorial term reduces to the configurational entropy introduced in (97),