Two-species hardcore reversible cellular automaton: matrix ansatz for dynamics and nonequilibrium stationary state

In this paper we study the statistical properties of a reversible cellular automaton in two out-of-equilibrium settings. In the first part we consider two instances of the initial value problem, corresponding to the inhomogeneous quench and the local quench. Our main result is an exact matrix product expression of the time evolution of the probability distribution, which we use to determine the time evolution of the density profiles analytically. In the second part we study the model on a finite lattice coupled with stochastic boundaries. Once again we derive an exact matrix product expression of the stationary distribution, as well as the particle current and density profiles in the stationary state. The exact expressions reveal the existence of different phases with either ballistic or diffusive transport depending on the boundary parameters.


Introduction
The determination of the macroscopic behavior of a system from its microscopic dynamics is a major challenge of statistical mechanics. One would like, for instance, to derive the hydrodynamic equations which govern the flow of conserved physical quantities (energy, electric charge, mass, etc.) starting from the knowledge of the microscopic interactions between the elementary constituents of the system. The problem is in general mathematically very difficult to address, because of the huge number of degrees of freedom of the many-body systems. In the last years the study of integrable systems and the development of the generalized hydrodynamics framework lead to the posibility of a precise large scale description of the dynamics of the charge and current profiles of all conserved quantities [1,2]. Nonetheless, the results of the generalized hydrodynamics rely on uncontrolled assumptions which can typically be corroborated only by the numerical checks. It is therefore important to investigate other theoretical techniques, which will enable us to study the non-equilibrium behavior of many-body systems analyticaly. Here we develop the matrix ansatz method to compute the time evolution of the probability distributions and to obtain the stationary states. The matrix ansatz technique has been first introduced in the context of Totally Asymmetric Simple Exclusion Process, to express the steady state [3], and has attracted a lot of attention since then in the context of continuous-time [4][5][6][7][8][9][10][11] and discretetime [12][13][14] classical stochastic processes and also in the context of open quantum systems [15,16]. In order to further develop the method, it thus appears primordial to design very simple models in which it is possible to perform analytical computations. One of the simplest class of many-body dynamical systems one can think of are the cellular automata. These are deterministic systems defined on a discrete space-time lattice with a discrete (typically finite) set of states at each space-time point. In particular, reversible cellular automata can be considered as minimal models for Hamiltonian (conservative) dynamics in the context of statistical mechanics [17,18]. In the present paper we investigate exact matrix product solutions for a specific 2-species reversible cellular automata introduced in [19,20] which may be considered as a model of a charged hardcore gas in one dimension. From the physical perspective the model is interesting since it exhibits a crossover between ballistic and diffusive transport.
Here we study the model in three different settings. The first one describes an inhomogeneous (bipartite) quench obtained by joining together two equilibrium distributions at different charge densities. The second one describes a local defect quench obtained by perturbing an equilibrium distribution on a single site. This setup can also be understood as a time propagation of the local observable. Less general quench setups were studied already in [19,20] using a different method. The last setting describes a boundary driving obtained by stochastic interactions of the system with the particle reservoirs.
The objective of the paper is twofold. Firstly, from a more mathematical perspective, we provide an exact matrix product expressions of the time-evolution of probability distributions using a unified framework, which relies on the algebraic "cancellation mechanism". It is one of the first examples of the explicit time-dependent matrix product ansatz beyond the recent results [21,22], which rely on the soliton counting. We expect that our framework can be generalized to more complicated systems. We also construct the exact matrix product expression of the stationary state of a stochastic boundary driven system. Secondly, from a more physical perspective, we provide new results associated with a phase transition between ballistic and diffusive phases of a boundary driven model.
In section 2 we present the dynamics of the model and introduce the different physical settings that we are interested in, namely the two quench protocols and the stochastic boundary driving. Then in section 3 (and in section 4) we study in details the inhomogeneous quench protocol (and the local quench protocol respectively). We provide an exact matrix product expression of the time-evolution of the system and we use it to compute analytically the density profiles. Finally, in section 5 we provide an exact matrix product expression of the stationary state of the boundary driven system. We compute analytically the particle currents and the density profiles for any size of the system. This allows us to identify the transition between the ballistic phase and the diffusive phase and obtain a complete phase diagram of the boundary driven system.

A two-component reversible cellular automaton 2.1 Definition of the models
We are interested in different models that share the same bulk dynamics, but with different boundary conditions and/or different initial conditions. In the first two cases, we will consider a periodic lattice with two different initial conditions (inhomogeneous and local quenches), while in the last case, we will be on a finite segment connected to two reservoirs. In the first setting, we will be interested in the propagation cone emerging from the origin, with the assumption that the propagation cone is not affected by the periodicity (i.e. we will always consider the case where the time t < L/4).
In order to give a precise mathematical definition of the models and to study them, we need to introduce some notations. The models are defined on a finite lattice comprising L sites labeled by an integer i ∈ [1, L]. Each site of the lattice is either empty or can carry a particle of species (or charge) + or −. We denote by τ t i the local occupation variable at site i and at time t. More precisely τ t i = 0 (respectively τ t i = +, τ t i = −) if there is a vacancy (respectively a particle of charge +, a particle of charge −) on site i at time t. The configuration of the lattice at time t is thus given by the L-uplet τ t = (τ t 1 , τ t 2 , . . . , τ t L ). Our aim is to study the statistical properties of the models. For this purpose it is usefull to introduce probability distributions. For τ 1 , . . . , τ L ∈ {0, +, −}, we denote by p t τ 1 ,...,τ L the probability that τ t = (τ 1 , τ 2 , . . . , τ L ). It will be convenient to encompass the probabilities of all configurations in a single vector where are the elementary basis vectors of C 3 . We now present the dynamics of the models. As already mentioned, the models that we are interested in share the same bulk dynamics. It is deterministic and is defined locally by updating a pair of neighboring sites (τ t+1 where τ ∈ {0, +, −} and , ∈ {+, −}. In words, the particles are propagating freely through the vacancies and there is a hard-core interaction between the particles, so that they cannot exchange their positions. We are considering a two-step discrete-time dynamics. During the first time step only the pairs of odd-even sites are updated whereas during the second time step only the pairs of even-odd sites are updated. This can be summarized as (τ t+1 for t − i even. For the periodic lattice, the complete dynamics is fully defined by the previous rules and the fact that the indices should be considered modulo the size of the system L. Note that in this case the number of sites should be even. For the open case, we are considering a lattice with an odd number of sites, so that at each time-step either the first or the last site of the lattice is singled out. When t is even we update τ t 1 with the following stochastic rule: τ t+1 1 is equal to 0 with probability α 0 , to + with probability α + , and to − with probability α − (where α 0 + α + + α − = 1). All other sites are updates accordingly to the bulk dynamics. Similarly in the case of odd t, we update τ t L with the following stochastic rule: τ t+1 L is equal to 0 with probability β 0 , to + with probability β + , and to − with probability β − (where β 0 + β + + β − = 1), and the remaining sites accordingly to the bulk dynamics.
In both periodic and open lattice settings the time evolution of the probability vector p(t) obeys a master equation of the following form where the two operators U e and U o correspond to the two different time-steps. In the first setting with periodic boundary conditions, we have The subscripts indicate on which sites of the lattice the operators are acting non-trivially. The operator U is a 9 × 9 matrix acting on two tensor space components and encoding the local update rule on two neighboring sites. It reads Note that it satisfies the braid relation leading to an R-matrix which obeys the braided Yang-Baxter equation. The fact that R(0) = U implies the integrability of the model with periodic boundary conditions. In the second, open lattice setting with stochastic boundary driving, we have Once again the subscripts indicate on which sites of the lattice the operators are acting nontrivially. The boundary operators B and B are 3 × 3 matrices acting on a single tensor space component and encoding the stochastic rules on the first and last sites. They read Note that they can not be obtained from solutions to the reflection equation associated with the R-matrix given in (8).

Separable states
We present here the building blocks of the quench protocols studied in sections 3 and 4. The bulk system dynamics allows for the existence of a class of spatially homogeneous currentcarrying states invariant under the system dynamics. They are defined for a system with periodic boundary conditions (the system size L is then even), or alternatively for a system with (fine-tuned) open boundary conditions (L is then odd). In both cases the respective stationary probability distribution obtains a factorized form, For this reason we shall call them separable states. It is easy to see that U o p = p , U e p = p, if and only if which amounts to the condition In case of the boundary driven system, with boundary matrices corresponding to (10), the distribution (11), (12) is stationary provided that the relation (14) holds. Solution (11), (12) then corresponds to the staggered profile of particle density, which fits the respective boundaries perfectly. For definiteness we shall parametrise β ± = µα ± , automatically satisfying (14). We also assume normalization α + + α − + α 0 = 1, β + + β − + β 0 = 1 throughout this section. Separable states are current carrying states. Indeed, the current of charges, defined as the average number of particles of a given charge crossing a bond, is (note that the factor 1/2 appears due to the fact that during one unit of time only one half of the bonds are involved in the dynamics). Below we shall describe several shock types involving the separable states. This will provide an intuitive physical understanding of the exact computations performed in sections 3, 4 and 5.
Shock between current-carrying separable state and separable state with no current A separable state satisfying (14), characterized by parameters α = (α 0 , α + , α − ) and β = (β 0 , β + , β − ), will be called p [i,j] (α, β), where [i, j] is the support of the state (i.e. the sublattice (interval) of sites on which the state is specified/localized). A special case of a separable states is the state p [i,j] (α, α) with µ = 1, i.e. α τ = β τ , τ = 0, ±. It corresponds to a constant density profile. A domain wall connecting p [−∞,i] (β, α) to p [i+1,∞] (β, β) has very simple dynamics: it remains sharp and it moves with velocity v = +1 or v = −1 depending on whether it is initially located at an odd or even site (due to the staggered nature of the propagator). Indeed, the action of the evolution operator on the respective domain wall is Likewise, the domain wall between p [−∞,i] (α, α) and p [i+1,∞] (β, α), is moving with the velocity ±1, depending on the initial position, Shock between two current-carrying separable states Let us join two separable states of type p [−∞,i] (α, γ), with γ ± = µα ± , on the left and p [i+1,∞] (δ, β), with δ ± = νβ ± , on the right. They carry respectively the currents Due to the current mismatch, an interface must be formed, the motion of which is regulated by the mass conservation condition. Supposing that the interface will form a stable shock, we find the interface velocity v + sh between + particles from the Rankine-Hugoniot condition [23] v + sh = Likewise, the interface velocity v − sh between − particles is Note that we used On the physical grounds we must have giving a consistency condition on the parameters of the separable states. Other consistency conditions are obtained requiring stability of the shock, according to a standard Lax theory of hydrodynamic shocks [23]. It is formulated in terms of characteristic velocities of infinitesimal perturbations on the top of a homogeneous background. The characteristic velocities are eigenvalues of the flux Jacobian DJ , = ∂J /∂ρ , for , = ±. Using (15), (21), for the left separable state we find J L ± = ρ L ± 1−µ 1+µ , and consequently the flux Jacobian is diagonal The corresponding characteristic velocities are equal c L 1 = c L 2 = 1−µ 1+µ . Similarly, for the right separable state we find DJ R , = ν−1 1+ν δ , , so that the characteristic velocities are

Out of equilibrium problems: quenches and boundary driven stationary state
In this subsection we present three non-equilibrium problems that we will study analytically. The first two are related to the study of inhomogeneous quench and local quench protocols. The third problem concerns the stationary state of the boundary driven model. In the case of the inhomogeneous quench, we study the situation where at t = 0 the system has a density of vacancies α 0 , the density of positively charged particles α + , and the density of negatively charged particles α − on the non-positive sites and densities β 0 , β + and β − on positive sites. The parameters satisfy α 0 +α + +α − = 1 and β 0 +β + +β − = 1. In mathematical terms the initial state (probability vector) is given by where the sites labelling should be understood modulo L (we have used negative labelling for later convenience). The initial state is obtained by joining two different stationary measures. We have indeed: implying that at time t only the sites −t < i ≤ t are non-trivially affected by the dynamics (we have a propagation cone) and the state of the system can be written as 1 In section 3 we provide an exact expression of the vector ψ(t) ∈ C 3 ⊗2t . It encodes nontrivial particle currents between the left and the right parts of the system. We use the exact expression to compute analytically the time-dependent density profiles. For the local quench, we are interested in the situation where at t = 0 the system has a density ρ 0 of vacancies, ρ + of positively charged particles and ρ − of negatively charged particles on all sites of the lattice except on site 1 where the densities are λ 0 , λ + and λ − , with the parameters satisfying the normalization conditions ρ 0 +ρ + +ρ − = 1 and λ 0 +λ + +λ − = 1. In this case the initial state (probability vector) takes the following form The initial state can be seen as a local defect at site 1 of an equilibrium state (i.e. of an invariant measure), since again the stationarity condition is satisfied. It implies that at time t only the sites −t < i ≤ t are non-trivially affected by the dynamics (we have a propagation cone) and the state of the system can be written In section 4 we provide an exact expression of the vector ψ(t) ∈ C 3 ⊗2t , which we then use to calculate the time-dependent density profiles.
The last out-of-equilibrium problem that we will consider is the computation of the stationary state of the model with stochastic boundary conditions that we introduced in the previous subsection. The operator U = U e U o indeed defines an irreducible Markov process that admits a unique stationary state p ∈ C 3 L , i.e. a unique probability vector satisfying Up = p. The probability vector p(t) is converging in the long-time limit toward this stationary state whatever the initial condition p(0). More precisely, we are looking for a pair of vector p, p satisfying This would obviously imply that p is the stationary state. In the section 5 we provide an exact construction of the vectors p and p . It allows us to compute analytically the particle currents and the density profiles in the stationary state.

Matrix product ansatz
Here we introduce the method used to construct the exact solutions to the three different out-of-equilibrium problems defined in the above subsection. The technique is called matrix ansatz and has been widely used in the statistical physics community to compute stationary states of stochastic boundary driven models. Here we will also use it to express exactly the time evolution of the probability distribution in the quench protocols. The idea is to express the components of the vector defined in (26) as a product of matrices (contracted with the row and the column vectors on the left and on the right to get a number) Two triples of matrices V τ , W τ , τ = 0, ±, are acting in an auxiliary space (which is different from the physical space of configurations), l| is a row vector of the auxiliary space and |r is a column vector of the auxiliary space. We will also use a similar matrix ansatz to construct exactly the solution of the local quench protocol defined in (29) and of the stationary state of the boundary driven model. The matrices V τ , W τ and the boundary vectors l| and |r have to be carefully chosen in order for this ansatz to correctly encode the probability distribution.
In the next sections we will provide an explicit expression of the matrices and of the boundary vectors and we will present an algebraic "cancellation scheme" to prove efficiently that the matrix ansatz is indeed correct.
We will now introduce all of the objects that will be needed to construct the matrix product expressions of the probability distributions studied in this paper. First of all we need to specify the auxiliary space on which the matrices and their building blocks will act and in which the boundary vectors will live. We define the Fock space F = span{|k } k≥0 and the creation and annihilation operators a, a † ∈ End(F), which satisfy the relation aa † = 1, where 1 is an identity operator on F. We also introduce the projection operator s = 1 − a † a = |0 0|, which satisfies the relations sa † = 0 and as = 0. It proves useful to introduce the coherent states where x is an arbitrary complex number. They satisfy Let us stress the difference between the state |k , k ∈ Z + , an element of the canonical basis (called from now on canonical state), and the coherent state |k , k ∈ C. They coincide only for k = 0. Note that all of these operators and vectors have already proven to be relevant in the context of exact matrix product expression of out-of-equilibrium stationary distributions in stochastic systems [3,5]. We also introduce a finite dimensional auxiliary vector space C 2 , with the canonical basis Additionally, it proves useful to introduce matrices A, B ∈ End(C 2 ) obeying BA = 0, A 2 = A and B 2 = B, which can be represented as the following 2 × 2 matrices: All these algebraic objects will be used extensively in the construction of the exact matrix product states in the following sections.

Inhomogeneous quench
This section is devoted to the study of the inhomogeneous quench protocol introduced above. We provide an exact construction of the time-dependent probability distribution in a matrix product form and we give an analytical expression of the density profiles. The exact results obtained prove rigorously the existence of both ballistic and diffusive transport.

Time-dependent matrix product ansatz
The vector ψ(t) defined in (26) encodes the non-trivial correlations induced by the inhomogeneous quench inside of the propagation cone. The main result of the present subsection is to provide an exact expression of its components (defined in equation (31)) using a matrix product ansatz The matrices V τ , W τ , τ = 0, ±, are operators acting on the auxiliary space C 2 ⊗ F, while the boundary vectors l| ∈ C 2 ⊗ F * and |r ∈ C 2 ⊗ F single out an element of contracted matrices which yields correct probability amplitudes ψ t τ −t+1 ,...,τt . For the matrices ( = ±), we have the following explicit matrix representation: and for the boundary vectors we have where we have used the canonical basis vectors e j , j = 1, 2 for C 2 and the B- We now show that the matrix product expression (38) indeed provides the correct timedependent probability distribution. The proof relies on a simple "telescopic scheme" based on three sets of algrebraic relations. In order to present efficiently those relations and the proof, we have to reformulate the matrix product expression using the tensor space formalism. The vector ψ(t) can be expressed as follows where V and W are the 3-component vectors encompassing the matrices The subscripts in (41) denote the components of the tensor space (i.e. the sites of the lattice) on which the vectors are located. Note that the components of the expression (41) are equivalent to (38). The first set of algebraic relations concerns the commutation relations between the matrices and can be concisely written using the tensor space formalism that is to say in components The other two sets of algebraic relations concern the boundary vectors on the left and on the right and read or equivalently These relations can be checked by direct computation.
The algebraic relations ensure that the propagation equation is fulfilled The matrix product structure of the time-dependent probability distribution could appear complicated at first sight but it turns out to be very efficient to compute analytically physical quantities such as density profiles.

Exact expression of physical observables
The first step toward the exact computation of physical observables is to verify that the matrix product construction of the above subsection defines a well-normalized probability distribution. The normalisation Z t is defined as the sum of the entries of the vector ψ(t). It can be checked from the explicit expression of the matrices (39) and boundary vectors (40) that ψ(t) is correctly normalized, i.e. where In order to lighten the notations, we have introduced the parameters α = α + + α − and β = β + + β − , which will be used throughout the paper. The details of the computation are given in appendix A.1.
The probability to observe a particle of species at site j = 2k − t, with 1 ≤ k ≤ t, can be easily expressed within the matrix product framework as This formula can be exactly evaluated using the explicit expression of the matrices and boundary vectors, see appendix A.2 (50) Note that the total particle density is constant and does not depend on time or on the position (as long as t − j is even). This is due to the fact that the total particle density is equal to one minus the density of vacancies. Indeed, the vacancies are propagating at velocity one (they do not interact) and hence their density at any time is trivially deduced from their density at initial time and is equal to 1 − α for t − j even.
The probability to observe a particle of species at site j = 2k + 1 − t, with 0 ≤ k ≤ t − 1, can be exactly computed using the matrix product expression, see appendix A.2 (52) Again, the total particle density has a very simple expression which does not depend on the time and on the position (as long as t − j is odd). Note that outside the propagation cone, for j ≤ −t or j > t, it is straightforward to compute

Hydrodynamic limit
This subsection is devoted to the study of the large time limit of the density profile. For this purpose we will use an equivalent expression of the density, see appendix A.2 for j = 2k − t, i.e. in the case where t − j is even. A similar expression also exists for The asymptotic behavior of the contour integral on the Euler scale j = xt, with −1 ≤ x ≤ 1, and t goes to infinity, can be studied using the saddle point analysis, which yields lim t→∞ n e xt, t = At the junction of the two density profiles x 0 = α−β α+β , where we observe a discontinuity on the Euler scale, the profile exhibits a smooth transition on the diffusive scale x−x 0 t 2 [19,20]. More precisely, for j = where we recall the definition of the error function Similarly, we have .
The equations (59) and (61) are proven starting from the contour integral expressions (55) and (56), using the expansion (β − αz) −1 = β −1 ∞ l=0 (zα/β) l . A saddle point analysis reveals that the terms of the sum with non-vanishing contribution in the limit t → ∞ are those for which l is of the order √ t, and yields the desired formulas. Note that the time evolution of the domain wall from the exact Matrix Product Ansatz (Fig. 3) shows a domain wall between the two separable states defined in subsection 2.2: p [−∞,i] (α, γ) with γ ± = β α α ± on the left and p [i+1,∞] (δ, β) with δ ± = α β β ± on the right. It is instructive to characterize it from the hydrodynamic point of view. From (19) and (22) in accordance with exact solution, and in addition c L 1 = c R 1 = v sh = α−β α+β . According to the Lax criterium, the shock in Fig.3 is marginally stable. It would be interesting to investigate other scenarios testing the hydrodynamic predictions, using Matrix Product Ansatz or direct Cellular Automata simulations.
We also observe that the asymptotic value of the density outside the propagation cone, i.e for x < −1 or x > 1, is given by It yields to discontinuities at the borders of the propagation cone moving with velocities −1 and 1 respectively. Those can also be interpreted as a domain walls between two separable states defined in subsection 2.2. The domain wall at the ray x = −1 is stable and characterized by the separable states p [−∞,i] (α, α) on the left and p [i+1,∞] (α, γ) with γ ± = β α α ± on the right. The domain wall at the ray x = 1 is stable and characterized by the separable states p [−∞,i] (β, δ) with δ ± = α β β ± on the left and p [i+1,∞] (β, β) on the right.

Local quench
This section is devoted to the study of the local quench protocol introduced above. Once again we provide an exact construction of the time-dependent probability distribution in a matrix product form and we give an analytical expression of the density profiles.

Time-dependent matrix product ansatz
The vector ψ(t) defined in (29) encodes the non-trivial correlations induced by the local quench in the propagation cone. The main result of the present subsection is to provide an exact expression of its components (defined in equation (31)) using a matrix product ansatz Note that the ansatz is a bit more complicated than previously because it involves a boundary vector l| τ which depends on the content of the site −t + 1. This dependence is reminiscent to the defect in the initial condition and could be equivalently encoded using a defect matrix (i.e different from V τ , W τ ) on the site −t + 1. The matrices V τ , W τ , τ = 0, ±, are again operators acting in the auxiliary space C 2 ⊗ F. The boundary vectors l| τ ∈ C 2 ⊗ F * and |r ∈ C 2 ⊗ F are used to contract the matrix product ansatz, thus yielding the components of the probability distribution. We choose the following explicit expression for the matrices Note that P is a projector (P 2 = P) obeying V τ P = V τ and W τ P = W τ , τ = ± , 0. For the boundary vectors, we choose the following form where we have used a coherent state in l| 0 , and canonical states in l| ± and |r . Note that l| τ P = l| τ , τ = ± , 0. In order to efficiently show that the matrix product expression (64) provides the correct time-dependent probability distribution, we reformulate the matrix product expression using the tensor space formalism where V , W and l| are the 3-component vectors The subscripts in (67) denote the components of the tensor space (i.e. the sites of the lattice) on which the vectors are located. The matrices obey the bulk relations (43). On the left boundary, the following condition should be satisfied that is to say On the right boundary, we have and finally, the boundary vectors should also satisfy the 'initial condition' The algebraic relations ensure that the propagation equation is fulfilled Once again the matrix product structure will prove to be very efficient to analytically compute physical quantities, such as the density profiles.

Exact expression of physical observables
The first step toward the exact computation of the time-dependent density profiles is to check that the vector ψ(t), defining the probability distribution in the propagation cone, is correctly normalized. The normalization Z t of the vector ψ(t) is defined as the sum of its components and can be written using the matrix product formalism as follows where It can be evaluated, see appendix B.1, using the explicit form of the matrices and of the boundary vectors where we have introduced in order to lighten the notations. Those parameters will also be used in the rest of the section. The probability to observe a particle of species at site j = 2k − t, with 1 ≤ k ≤ t, can be easily expressed within the matrix product framework as This formula can be exactly evaluated using the explicit expression of the matrices and boundary vectors, see appendix B.2 The total particle density takes the very simple expression The probability to observe a particle of species at site j = 2k + 1 − t, with 1 ≤ k ≤ t − 1, can be exactly computed using the matrix product expression, see appendix B.2  Once again, the total particle density is simply The density on the site j = −t + 1 (i.e. for k = 0) takes a specific expression, see appendix B.2 In this case the total particle density reads

Hydrodynamic limit
An equivalent alternative expression for the density is given by Using the saddle point analysis we can again study the asymptotic behavior of the contour integral in the regime where j = xt, with −1 ≤ x ≤ 1, and t goes to infinity. The contribution of the contour integral is exponentially small when x = 0, and one gets lim t→∞ n e xt, t = lim t→∞ n o xt, t = ρ .
Note that at the site j = −t + 1, we can obtain an explicit expression, which in the long time limit reads lim This corresponds to the asymptotic value of a ballistically propagating peak, which can be observed in Fig. 4. In the vicinity of x = 0, in the regime x = y √ t , we obtain The asymptotic behavior of n o y √ t, t is exactly the same. The ballistic left mover on Fig. 4 on the extremity of the propagation cone can be interpreted as a consequence of the ballistically propagating vacancies. At time t = 0 there is a discrepancy between the density at site 1 and the density elsewhere. The latter is propagating ballistically to the left because of the staggered nature of the dynamics (if the local defect was initially located on an even site the ballistic peak would be moving to the right).

Matrix product expression of the stationary state
We recall that we are using the short-hand notations The main result of the present subsection is an exact expression of the components of the stationary states p, defined by the equations Up = p (where the operator U has been defined in the subsection 2.3). We have where the matrices V τ , W τ ∈ End(C 2 ⊗ F ⊗ F). Note that the auxiliary space is now a bit more complicated than in the previous cases because it involves two copies of the infinite dimensional Fock space F. The matrix product ansatz is constructed by employing matrices A and B, which were introduced in the preceding section (37): The vectors l| τ ∈ C 2 ⊗ F * ⊗ F * read explicitly The vectors |rr τ,τ ∈ C 2 ⊗ F ⊗ F are expressed using the B-eigenvector v = e 1 + e 2 (93) In order to prove the matrix product expression of the stationary state it will be useful to reformulate it using the tensor space formalism where We recall that the subscripts indicate on which tensor space components (i.e. on which sites of the lattice) the operators are acting. As already mentioned, the proof relies on the construction of a pair of vectors p and p satisfying the following relations These relations imply the following eigenvector relation We already provided the expression of the stationary state p. We now give the expression of the vector p . For this purpose we need to introduce the boundary vectors where and ll| 0τ = β 0 α τ α e t 1 ⊗ 0|a ⊗ β/β 0 |s + 1 + Building on those, we define The relations (96) are proven using local exchange relations in the bulk on the right boundary and on the left boundary By direct computation it is possible to check that all of the above relations are satisfied by the ansatz (91), (92) and (93).

Exact expression of physical observables
In order to compute physical quantities we need to evaluate the normalization of the matrix product probability distribution, which corresponds to the sum over all the entries (i.e over all configurations) of the stationary matrix product state where We show in appendix C.1 that it is given by the following exact formula The mean particle current associated to the species = ± is defined by the mean number of particles of species that jump between sites i and i + 1 during a unit of time. Since the bulk dynamics conserves the particle number, the steady state current does not depend on the site index. We can, for instance, compute it between site 1 and site 2. Using the matrix product expression of the stationary state, the current is given by The above expression can be evaluated exactly (see appendix C.2) From this result we can deduce that the total particle current takes a very simple form When α = β, the current is given by the appropriate limit of the general expression The mean particle density associated to the species = ± at site i is defined by the probability to observe a particle at site i in the stationary state. Since we are interested in mean behavior, we need to average the density over two time steps (odd and even): the stationary state is described by p and p respectively. The density at site i = 2k + 1 is thus given in the matrix product formalism by whereZ L is the normalization (i.e. the sum of the components) of p . Once again the matrix product expression can be computed explicitly (see appendix C.3) We deduce that the total particle density takes a very simple form Note that we have an equivalent expression of the particle density (205) which allows us to easily take the limit β → α, and obtain an expression of the particle density in the case α = β Note that in this case the density profile is a linear profile interpolating between the densities of the left and right reservoirs α and β .

Hydrodynamic limit and phase transitions
In this subsection we investigate the properties of the physical quantities in the large system size limit L → ∞. Depending on the boundary parameters, we encounter different scaling behaviors of the particle current with respect to the system size, which are interpreted as signature of ballistic and diffusive transport. More precisely the parameters space splits into three different phases: The right reservoir phase, for α < β. The particle current has a finite large system size limit lim The limit of the particle density is It has a constant value 1 2 β β (α + β), fixed by the right reservoir in the bulk, and an exponential tail near the left boundary.
The left reservoir phase, for β < α. The particle current has a finite large system size limit lim The limit of the particle density is It has a constant value 1 2 α α (α + β), fixed by the left reservoir, and an exponential tail near the right boundary.
The diffusive phase, for α = β. In this phase the particle current decays as 1/L in the large system size limit. More precisely we have The limit of the particle density is In this phase the Fick law (which is a signature of the diffusive behavior) is satisfied where D = 1−α α is the diffusion constant, in accordance with the value computed in [19,20].

Conclusion
In the article we presented a unified algebraic framework for obtaining exact solutions of different out-of-equilibrium setups in terms of the (time-dependent) matrix product ansatz.
In particular, using the new framework we obtained an explicit solution of the hardcore interacting deterministic lattice gas, and calculated physically interesting quantities such as the current and the density profiles. Interestingly, in the case of the boundary driving we observed that, depending on the parameters of external driving, the system undergoes a phase transition. Our framework paves a new way for obtaining exact time-dependent results in statistical systems (in particular, with deterministic bulk dynamics). We hope that the approach presented here can be generalized to more complicated systems. The first step on this path would be to understand the time-dependent solution of the classical and quantum Rule 54 automaton [21,22].
The full time dependent result also enables the calculation of physical quantities beyond the expectation values of currents or charges. In particular it would be interesting to obtain a large deviation functional, especially in the case of the boundary driving, where we observed the phase transition.

A.2 Computation of the densities
The probability to observe a particle of species at site j = 2k − t, with 1 ≤ k ≤ t, is given by the matrix product expression Using the commutation property (124) and the eigenvectors relations (126) we can simplify the expression The next step is to use the closure relation in the Fock space and insert it to the left of W . Note that the integration contour includes the pole at z = 0 and excludes all others. Then using the fact that the coherent states |z and 1/z| are respectively eigenstates of the annihilation and creation operators a and a † , we obtain (131) The integrand can be evaluated using the triangular structure of matrices A and B The residue at z = 0 corresponding to the first of the two previous terms is equal to whereas the residue at z = 0 corresponding to the second term is equal to Gathering these results together we end up with the following expression n e (j, t) = β β Finally, using the fact that α 0 + α = 1 and β 0 + β = 1, the expression reduces to (136) Note that we also have an alternative expression for the density, if we do not evaluate explicitly the residue (133) where the contour of the integral encircles the pole at z = 0 and excludes the pole at z = β/α. The probability to observe a particle of species at site j = 2k + 1 − t, with 0 ≤ k ≤ t − 1, is given by the matrix product expression A very similar computation yields Using the fact that α 0 + α = 1 and β 0 + β = 1, the expression reduces to (140) As before, we also have a contour integral expression B Details of the computations for the local quench

B.1 Computation of the normalisation
The normalisation of the time-dependent matrix product state is defined as where To evaluate this quantity explicitly we first use the fact that which can be readily checked using the explicit expression of the matrices (65). The normalization can then be rewritten as The boundary vectors l| τ and |r (see explicit expressions in (66)) fulfill the following eigenvalue relations We have also the scalar product expression Those relations lead immediately to the expression

B.2 Computation of the densities
The probability to observe a particle of species at site j = 2k −t, with 1 ≤ k ≤ t, is expressed within the matrix product framework as Using the commutation property (143) and the eigenvectors properties (145) the expression is simplified into Then using the explicit expression of the matrices (65) and of the boundary vectors (66) and the properties V τ P = V τ , W τ P = W τ and l| τ P = l| τ , τ = ± , 0, we can get rid of the projector P in the matrix product and we obtain (150) The next step is to use the closure relation in the Fock space and insert it once between 2 the operators a and a † . Note that the integration contour includes the pole at z = 0 and excludes all others. Then using the fact that the coherent states |z and 1/z| are respectively eigenstates of the annihilation and creation operators a and a † , we obtain 2 you can indeed observe in equation (150) that the operators a are all located on the left side and the operators a † are all located on the right side of the expression If we do not evaluate explicitly the residue, we obtain an alternative expression of the densities, Finally we have the particular case of the density of species at site j = −t + 1. The corresponding matrix product expression is Using the commutation property (143), the right eigenvector property in (145) we obtain Then using recursively the relations we get The scalar product relations l| · |r = λ , yield the final expression Taking into account the fact that ρ 0 + ρ = 1 and λ 0 + λ = 1, the formula reduces to C Details of the computations for the boundary driven model
Taking into account that α + α 0 = 1 and β + β 0 = 1, the expression reduces to We now come back to the proof of equation (184). The first step is to observe that the action of the coherent state |α/α 0 in the second tensor component of the right boundary vectors can be evaluated, similarly as in the previous computations. The action of the vector 0| in the third tensor component of the left vector can be directly evaluated on (W 0 + W + + W − ) n . This last step kills the term involving the matrix B in W 0 + W + + W − (see equation (172)), so that e t 1 is a left eigenvector of the remaining terms (e t 1 is a left eigenvector of the matrix A). Applying all these manipulations we end up with