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

TL;DR

This paper analyzes a reversible two-species cellular automaton's dynamics and stationary states using an exact matrix product approach, revealing phase-dependent transport behaviors.

Contribution

It introduces an exact matrix product solution for the automaton's dynamics and stationary states, providing analytical insights into phase-dependent transport phenomena.

Findings

Exact matrix product expressions for time evolution and stationary distribution

Identification of phases with ballistic or diffusive transport

Analytical density profiles and particle currents

Abstract

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.

Figures

Equations

$$\mathbf{p}(t)=\sum _{{\tau }_i\in \{0,+,-\}}{p}_{{\tau }_1,\dots ,{\tau }_L}^t{e}_{{\tau }_1}\otimes e_{{\tau }_2}\otimes \cdots \otimes e_{{\tau }_L}$$

$$e_0=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),e_{+}=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),e_{-}=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$$

$$\varphi (0,\tau )=(\tau ,0),\varphi (\tau ,0)=(0,\tau ),\varphi (ϵ,{ϵ}^{\prime })=(ϵ,{ϵ}^{\prime })$$

$$\mathbf{p}(t+1)=\{\begin{array}{l}{\mathbb{U}}^{\mathsf{e}}\mathbf{p}(t),t\text{mbox}even,\\ {\mathbb{U}}^{\mathsf{o}}\mathbf{p}(t),t\text{mbox}odd,\end{array}$$

$${\mathbb{U}}^{\mathsf{e}}=\prod _{k=1}^{\frac{L}{2}}{U}_{2k,2k+1}\text{mbox}and{\mathbb{U}}^{\mathsf{o}}=\prod _{k=1}^{\frac{L}{2}}{U}_{2k-1,2k}.$$

$$U=\left(\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right).$$

$$U_{12}{U}_{23}{U}_{12}=U_{23}{U}_{12}{U}_{23},$$

$$R(z)=P\frac{\mathbb{I}+zU}{z+1}$$

$$t_{\lambda}(z)=tr_{0}\Big{(}R_{0L}(z)\,R_{0,L-1}(z+\lambda)\,R_{0,L-2}(z)\,R_{0,L-3}(z+\lambda)\,\cdots R_{02}(z)\,R_{01}(z+\lambda)\Big{)},$$

$$t_{\lambda }(-\lambda {)}^{-1}{t}_{\lambda }(0)=V_{12}^{+}{V}_{34}^{+}\cdots V_{L-1,L}^{+}{V}_{L-2,L-1}^{+}{V}_{L-4,L-3}^{+}\cdots V_{23}^{+}{V}_{L,1}^{+},$$

$$\lim_{\lambda\to\infty}\big{(}t_{\lambda}(-\lambda)^{-1}\,t_{\lambda}(0)\big{)}={\mathbb{U}}^{\mathsf{o}}{\mathbb{U}}^{\mathsf{e}}$$

$$t_{\infty}(z)=tr_{0}\Big{(}R_{0L}(z)\,\check{U}_{0,L-1}\,R_{0,L-2}(z)\,\check{U}_{0,L-3}\,\cdots R_{02}(z)\,\check{U}_{01}\Big{)},\quad\mbox{where}\quad\check{U}=PU.$$

$${\mathbb{U}}^{\mathsf{o}}=\left(\prod _{k=1}^{\frac{L-1}{2}}{U}_{2k-1,2k}\right){\overline{B}}_L\text{mbox}and{\mathbb{U}}^{\mathsf{e}}=B_1\left(\prod _{k=1}^{\frac{L-1}{2}}{U}_{2k,2k+1}\right).$$

$$B=\left(\begin{array}{ccc}{\alpha }_0& {\alpha }_0& {\alpha }_0\\ {\alpha }_{+}& {\alpha }_{+}& {\alpha }_{+}\\ {\alpha }_{-}& {\alpha }_{-}& {\alpha }_{-}\end{array}\right),\overline{B}=\left(\begin{array}{ccc}{\beta }_0& {\beta }_0& {\beta }_0\\ {\beta }_{+}& {\beta }_{+}& {\beta }_{+}\\ {\beta }_{-}& {\beta }_{-}& {\beta }_{-}\end{array}\right).$$

$$p_{s_1{s}_2{s}_3\dots }$$

$$p_{s_1{s}_2{s}_3\dots }^{\prime }$$

$$U\left(\begin{array}{c}{\alpha }_0\\ {\alpha }_{+}\\ {\alpha }_{-}\end{array}\right)\otimes \left(\begin{array}{c}{\beta }_0\\ {\beta }_{+}\\ {\beta }_{-}\end{array}\right)=\left(\begin{array}{c}{\beta }_0\\ {\beta }_{+}\\ {\beta }_{-}\end{array}\right)\otimes \left(\begin{array}{c}{\alpha }_0\\ {\alpha }_{+}\\ {\alpha }_{-}\end{array}\right),$$

$$\frac{{\alpha }_{+}}{{\alpha }_{-}}=\frac{{\beta }_{+}}{{\beta }_{-}}.$$

$$J_{+}$$

$$J_{-}$$

$$\mathbb{U}{\mathbf{p}}_{[-\infty ,0]}(\mathbf{\beta },\mathbf{\alpha })\otimes {\mathbf{p}}_{[1,\infty ]}(\mathbf{\beta },\mathbf{\beta })$$

$$\mathbb{U}{\mathbf{p}}_{[-\infty ,-1]}(\mathbf{\beta },\mathbf{\alpha })\otimes {\mathbf{p}}_{[0,\infty ]}(\mathbf{\beta },\mathbf{\beta })$$

$$\mathbb{U}{\mathbf{p}}_{[-\infty ,0]}(\mathbf{\alpha },\mathbf{\alpha })\otimes {\mathbf{p}}_{[1,\infty ]}(\mathbf{\beta },\mathbf{\alpha })$$

$$\mathbb{U}{\mathbf{p}}_{[-\infty ,1]}(\mathbf{\alpha },\mathbf{\alpha })\otimes {\mathbf{p}}_{[2,\infty ]}(\mathbf{\beta },\mathbf{\alpha })$$

$$J_{+}^L$$

$$J_{+}^R$$

$$v_{sh}^{+}=\frac{J_{+}^R-J_{+}^L}{{\rho }_{+}^R-{\rho }_{+}^L}=\frac{{\beta }_{+}(\nu -1)-{\alpha }_{+}(1-\mu )}{{\beta }_{+}(1+\nu )-{\alpha }_{+}(1+\mu )}.$$

$$v_{sh}^{-}=\frac{J_{-}^R-J_{-}^L}{{\rho }_{-}^R-{\rho }_{-}^L}=\frac{{\beta }_{-}(\nu -1)-{\alpha }_{-}(1-\mu )}{{\beta }_{-}(1+\nu )-{\alpha }_{-}(1+\mu )}.$$

$${\rho }_{\pm }^L$$

$${\rho }_{\pm }^R$$

$$v_{sh}^{-}=v_{sh}^{+}=v_{sh},$$

$$\frac{\nu -1}{1+\nu }\le v_{sh}\le \frac{1-\mu }{1+\mu }.$$

$$\mathbf{p}(0)=\underset{-L/2<i\le 0}{⨂}\left(\begin{array}{c}{\alpha }_0\\ {\alpha }_{+}\\ {\alpha }_{-}\end{array}\right)\otimes \underset{0<i\le L/2}{⨂}\left(\begin{array}{c}{\beta }_0\\ {\beta }_{+}\\ {\beta }_{-}\end{array}\right),$$

Full text

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

M. Medenjak1, V. Popkov2,3,4, T. Prosen2, E. Ragoucy5, M. Vanicat2*

1 Institut de Physique Théorique Philippe Meyer, Ecole Normale Supérieure, PSL University, Sorbonne Universités, CNRS, 75005 Paris, France

2 Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia

3 Institut für Theoretische Physik, Universität zu Köln Zülpicher Str. 77, 50937 Köln

4 Bergisches Universität Wuppertal, Gauss Str. 20, D-42097 Wuppertal

5 Laboratoire d’Annecy-le-Vieux de Physique Théorique LAPTh, Université Grenoble Alpes, USMB, CNRS, F-74000 Annecy

• matthieu.vanicat@fmf.uni-lj.si

Abstract

**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. **

Contents

1. 1 Introduction

2. 2 A two-component reversible cellular automaton

3. 2.1 Definition of the models

4. 2.2 Separable states

5. 2.3 Out of equilibrium problems: quenches and boundary driven stationary state

6. 2.4 Matrix product ansatz

7. 3 Inhomogeneous quench

8. 3.1 Time-dependent matrix product ansatz

9. 3.2 Exact expression of physical observables

10. 3.3 Hydrodynamic limit

11. 4 Local quench

12. 4.1 Time-dependent matrix product ansatz

13. 4.2 Exact expression of physical observables

14. 4.3 Hydrodynamic limit

15. 5 Stochastic boundary driving

16. 5.1 Matrix product expression of the stationary state

17. 5.2 Exact expression of physical observables

18. 5.3 Hydrodynamic limit and phase transitions

19. 6 Discussion and conclusion

20. A Details of the computations for the inhomogeneous quench

21. A.1 Computation of the normalisation

22. A.2 Computation of the densities

23. B Details of the computations for the local quench

24. B.1 Computation of the normalisation

25. B.2 Computation of the densities

26. C Details of the computations for the boundary driven model

27. C.1 Computation of the normalisation

28. C.2 Computation of the particle currents

29. C.3 Computation of the densities

1 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 analytically. 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 the 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 discrete-time [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 a 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.

2 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\in[1,L]$. Each site of the lattice is either empty or can carry a particle of species (or charge) $+$ or $-$. We denote by $\tau_{i}^{t}$ the local occupation variable at site $i$ and at time $t$. More precisely $\tau_{i}^{t}=0$ (respectively $\tau_{i}^{t}=+$, $\tau_{i}^{t}=-$) 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 $\underline{\tau}^{t}=(\tau_{1}^{t},\tau_{2}^{t},\dots,\tau_{L}^{t})$.

Our aim is to study the statistical properties of the models. For this purpose it is useful to introduce probability distributions. For $\tau_{1},\dots,\tau_{L}\in\{0,+,-\}$, we denote by $p^{t}_{\tau_{1},\dots,\tau_{L}}$ the probability that $\underline{\tau}^{t}=(\tau_{1},\tau_{2},\dots,\tau_{L})$. It will be convenient to encompass the probabilities of all configurations in a single vector

[TABLE]

where

[TABLE]

are the elementary basis vectors of ${\mathbb{C}}^{3}$.

We now present the dynamics of the models. As already mentioned, the models that we are interested in share the same deterministic bulk dynamics. It is defined locally by updating a pair of neighboring sites $(\tau_{i}^{t+1},\tau_{i+1}^{t+1})=\phi(\tau_{i}^{t},\tau_{i+1}^{t})$ with the following rules

[TABLE]

where $\tau\in\{0,+,-\}$ and $\epsilon,\epsilon^{\prime}\in\{+,-\}$. 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 $(\tau_{i}^{t+1},\tau_{i+1}^{t+1})=\phi(\tau_{i}^{t},\tau_{i+1}^{t})$ 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 $\tau_{1}^{t}$ with the following stochastic rule: $\tau_{1}^{t+1}$ is equal to [math] with probability $\alpha_{0}$, to $+$ with probability $\alpha_{+}$, and to $-$ with probability $\alpha_{-}$ (where $\alpha_{0}+\alpha_{+}+\alpha_{-}=1$). All other sites are updated accordingly to the bulk dynamics. Similarly, if $t$ is odd, we update $\tau_{L}^{t}$ with the following stochastic rule: $\tau_{L}^{t+1}$ is equal to [math] with probability $\beta_{0}$, to $+$ with probability $\beta_{+}$, and to $-$ with probability $\beta_{-}$ (where $\beta_{0}+\beta_{+}+\beta_{-}=1$), and the remaining sites accordingly to the bulk dynamics.

In both periodic and open lattice settings the time evolution of the probability vector $\bm{p}(t)$ obeys a master equation of the following form

[TABLE]

where the two operators ${\mathbb{U}}^{\mathsf{e}}$ and ${\mathbb{U}}^{\mathsf{o}}$ correspond to the two different time-steps. In the first setting with periodic boundary conditions, we have

[TABLE]

The subscripts indicate on which sites of the lattice the operators are acting non-trivially. The operator $U$ is a $9\times 9$ matrix acting on two tensor space components and encoding the local update rule on two neighboring sites. It reads

[TABLE]

Note that it satisfies the braid relation

[TABLE]

leading to an $R$-matrix

[TABLE]

where $P$ is the permutation operator $P(u\otimes v)\equiv v\otimes u$ acting on the tensor space ${\mathbb{C}}^{3}\otimes{\mathbb{C}}^{3}$. The $R$-matrix obeys the Yang-Baxter equation (with additive spectral parameter $z$), the unitarity condition $R(z)R(-z)={\mathbb{I}}$ and $R(0)=PU$.

From this $R$-matrix and any $\lambda\in{\mathbb{C}}$, one builds the transfer matrix (we remind that $L$ is even for periodic boundary conditions)

[TABLE]

which commutes for different values of the spectral parameter. Then, standard calculation leads to

[TABLE]

where $V^{\pm}=\frac{{\mathbb{I}}\pm\lambda\,U}{1\pm\lambda}$. We have used the relations $(V^{-})^{-1}=V^{+}$ and $U^{2}={\mathbb{I}}$. From the property $\lim_{\lambda\to\infty}V^{\pm}=U$, one gets

[TABLE]

ensuring the integrability of the model with periodic boundary condition. ${\mathbb{U}}$ commutes with the transfer matrix

[TABLE]

In the second, open lattice setting with stochastic boundary driving, we have

[TABLE]

Once again the subscripts indicate on which sites of the lattice the operators are acting non-trivially. The boundary operators $B$ and $\overline{B}$ are $3\times 3$ matrices acting on a single tensor space component and encoding the stochastic rules on the first and last sites. They read

[TABLE]

Note that they can not be obtained from solutions to the reflection equation associated with the $R$-matrix given in (8).

2.2 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 current-carrying 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,

[TABLE]

where $\bm{p}$ and $\bm{p}^{\prime}$ represent probability vectors over all particle configurations after odd and even times, respectively. For this reason we shall call them separable states. It is easy to see that ${\mathbb{U}}^{\mathsf{o}}\bm{p}=\bm{p^{\prime}},{\mathbb{U}}^{\mathsf{e}}\bm{p^{\prime}}=\bm{p}$, if and only if

[TABLE]

which amounts to the condition

[TABLE]

In case of the boundary driven system, with boundary matrices corresponding to (13), the distribution (14), (15) is stationary provided that the relation (17) holds. Solution (14), (15) then corresponds to the staggered profile of particle density, which fits the respective boundaries perfectly. For definiteness we shall parametrise $\beta_{\pm}=\mu\alpha_{\pm}$, automatically satisfying (17). We also assume normalization $\alpha_{+}+\alpha_{-}+\alpha_{0}=1,\ \beta_{+}+\beta_{-}+\beta_{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

[TABLE]

where $\langle\tau\tau^{\prime}\rangle_{12}$ is the joint probability of having a $\tau$ particle at site 1 and $\tau^{\prime}$ particle at site 2, $\tau,\tau^{\prime}=0,+,-$. 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 (17), characterized by parameters $\bm{\alpha}=(\alpha_{0},\alpha_{+},\alpha_{-})$ and $\bm{\beta}=(\beta_{0},\beta_{+},\beta_{-})$, will be called $\bm{p}_{[i,j]}(\bm{\alpha},\bm{\beta})$, 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 $\bm{p}_{[i,j]}(\bm{\alpha},\bm{\alpha})$ with $\mu=1$, i.e. $\alpha_{\tau}=\beta_{\tau}$, $\tau=0,\pm$. It corresponds to a constant density profile. A domain wall connecting $\bm{p}_{[-\infty,i]}(\bm{\beta},\bm{\alpha})$ to $\bm{p}_{[i+1,\infty]}(\bm{\beta},\bm{\beta})$ 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

[TABLE]

We note here that in the context of hydrodynamics, such microscopically sharp moving interfaces are named contact discontinuities, and in our model they appear due to the free propagation of vacancies in the (disordered) environment of $+$ and $-$ particles with velocity $v=\pm 1$, i.e. the light-cone velocity. The respective normal modes are thus interaction-free modes, and their appearance can be seen in all quench scenarios considered in this paper, see Figs. 3,4, at the borders of the light cone $x=\pm t$.

Likewise, the domain wall between $\bm{p}_{[-\infty,i]}(\bm{\alpha},\bm{\alpha})$ and $\bm{p}_{[i+1,\infty]}(\bm{\beta},\bm{\alpha})$, is moving with the velocity $\pm 1$, depending on the initial position,

[TABLE]

Shock between two current-carrying separable states

Let us join two separable states of type $\bm{p}_{[-\infty,i]}(\bm{\alpha},\bm{\gamma})$, with $\gamma_{\pm}=\mu\alpha_{\pm}$, on the left and $\bm{p}_{[i+1,\infty]}(\bm{\delta},\bm{\beta})$, with $\delta_{\pm}=\nu\beta_{\pm}$, on the right. They carry the currents

[TABLE]

respectively. Due to the current mismatch, an interface or two separate interfaces must be formed, corresponding to shocks or rarefaction waves, or both (due to existence of two conservation laws for

• and - particles). Their motion is regulated by the mass conservation condition. Supposing that the interface will form a single shock between the two current carrying states, we find the interface velocity $v_{sh}^{+}$ between $+$ particles from the Rankine-Hugoniot condition [23]

[TABLE]

Likewise, the interface velocity $v_{sh}^{-}$ between $-$ particles is

[TABLE]

Note that we used

[TABLE]

For a single shock scenario to be consistent we must have

[TABLE]

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_{\epsilon,\epsilon^{\prime}}=\partial J_{\epsilon}/\partial\rho_{\epsilon^{\prime}}$, for $\epsilon,\epsilon^{\prime}=\pm$. Using (18), (24), for the left separable state we find $J^{L}_{\pm}=\rho^{L}_{\pm}\frac{1-\mu}{1+\mu}$, and consequently the flux Jacobian is diagonal $DJ_{\epsilon,\epsilon^{\prime}}^{L}=\frac{1-\mu}{1+\mu}\delta_{\epsilon,\epsilon^{\prime}}$. The corresponding characteristic velocities are equal $c_{1}^{L}=c_{2}^{L}=\frac{1-\mu}{1+\mu}$. Similarly, for the right separable state we find $DJ_{\epsilon,\epsilon^{\prime}}^{R}=\frac{\nu-1}{1+\nu}\delta_{\epsilon,\epsilon^{\prime}}$, so that the characteristic velocities are $c_{1}^{R}=c_{2}^{R}=\frac{\nu-1}{1+\nu}$. The Lax criterium for shock stability requires $\max(c_{1}^{L},c_{2}^{L})\geq v_{sh}$, $\min(c_{1}^{R},c_{2}^{R})\leq v_{sh}$, thus giving

[TABLE]

The condition (26) guarantees that small perturbations generated at both sides of a shock are absorbed by the shock interface, and therefore shock remains stable. Violating (26) results in a smoothening of the interface and to a formation of a rarefaction wave, thus leading to shock instability. The case when one or two characteristic velocities coincide with the shock velocity, is at the border between the two possible scenarios, and will be referred to as marginally stable shock situation.

2.3 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 $\alpha_{0}$, the density of positively charged particles $\alpha_{+}$, and the density of negatively charged particles $\alpha_{-}$ on the non-positive sites and densities $\beta_{0}$, $\beta_{+}$ and $\beta_{-}$ on positive sites. The parameters satisfy $\alpha_{0}+\alpha_{+}+\alpha_{-}=1$ and $\beta_{0}+\beta_{+}+\beta_{-}=1$. In mathematical terms the initial state (probability vector) is given by

[TABLE]

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. Indeed, we have:

[TABLE]

implying that at time $t$ only the sites $-t<i\leq t$ are non-trivially affected by the dynamics (we have a propagation cone) and the state of the system can be written as111For the sake of simplicity we omit to describe the second propagation cone emerging from site $L/2$. We recall that we will always consider the case where the time $t<L/4$, i.e. when the two propagation cones did not merge yet (the system then behaves as it was on an infinite line).

[TABLE]

In section 3 we provide an exact expression of the vector $\bm{\psi}(t)\in\left(\mathbb{C}^{3}\right)^{\otimes 2t}$. It encodes non-trivial 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 $\rho_{0}$ of vacancies, $\rho_{+}$ of positively charged particles and $\rho_{-}$ of negatively charged particles on all sites of the lattice except on site $1$ where the densities are $\lambda_{0}$, $\lambda_{+}$ and $\lambda_{-}$, with the parameters satisfying the normalization conditions $\rho_{0}+\rho_{+}+\rho_{-}=1$ and $\lambda_{0}+\lambda_{+}+\lambda_{-}=1$. In this case the initial state (probability vector) takes the following form

[TABLE]

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

[TABLE]

is satisfied. It implies that at time $t$ only the sites $-t<i\leq t$ are non-trivially affected by the dynamics (we have a propagation cone) and the state of the system can be written

[TABLE]

In section 4 we provide an exact expression of the vector $\bm{\psi}(t)\in\left(\mathbb{C}^{3}\right)^{\otimes 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 ${\mathbb{U}}={\mathbb{U}}^{\mathsf{e}}{\mathbb{U}}^{\mathsf{o}}$ indeed defines an irreducible Markov process that admits a unique stationary state $\bm{p}\in\left(\mathbb{C}^{3}\right)^{L}$, i.e. a unique probability vector satisfying ${\mathbb{U}}\bm{p}=\bm{p}$. The probability vector $\bm{p}(t)$ is converging in the long-time limit toward this stationary state for any initial condition $\bm{p}(0)$. More precisely, we are looking for a pair of vector $\bm{p},\bm{p^{\prime}}$ satisfying

[TABLE]

This would obviously imply that $\bm{p}$ is the stationary state. In the section 5 we provide an exact construction of the vectors $\bm{p}$ and $\bm{p^{\prime}}$. It allows us to analytically compute the particle currents and the density profiles in the stationary state.

2.4 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 (29)

[TABLE]

as a product of matrices (contracted with the row and the column vectors on the left and on the right to get a number)

[TABLE]

Two triples of matrices $V_{\tau}$, $W_{\tau}$, $\tau=0,\pm$, are acting in an auxiliary space (which is different from the physical space of configurations), ${\langle l|}$ is a row vector of the auxiliary space and ${|r\rangle}$ 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 (32) and of the stationary state of the boundary driven model. The matrices $V_{\tau}$, $W_{\tau}$ and the boundary vectors ${\langle l|}$ and ${|r\rangle}$ 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 ${\cal F}=\mbox{span}\{{|k\rangle}\}_{k\geq 0}$ and the creation and annihilation operators $\mathfrak{a},\mathfrak{a}^{\dagger}\in\mbox{End}({\cal F})$,

[TABLE]

which satisfy the relation $\mathfrak{a}\mathfrak{a}^{\dagger}=1$, where $1$ is an identity operator on ${\cal F}$. We also introduce the projection operator $\mathfrak{s}=1-\mathfrak{a}^{\dagger}\mathfrak{a}={|0\rangle}{\langle 0|}$, which satisfies the relations $\mathfrak{s}\mathfrak{a}^{\dagger}=0$ and $\mathfrak{a}\mathfrak{s}=0$. It proves useful to introduce the coherent states

[TABLE]

where $x$ is an arbitrary complex number. They satisfy

[TABLE]

Let us stress the difference between the state ${|k\rangle}$, $k\in{\mathbb{Z}}_{+}$, an element of the canonical basis (called from now on canonical state), and the coherent state ${|\underline{k}\rangle}$, $k\in{\mathbb{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 ${\mathbb{C}}^{2}$, with the canonical basis

[TABLE]

Additionally, it proves useful to introduce matrices $A,B\in{\rm End(}{\mathbb{C}}^{2})$ obeying $BA=0$, $A^{2}=A$ and $B^{2}=B$, which can be represented as the following $2\times 2$ matrices:

[TABLE]

All these algebraic objects will be used extensively in the construction of the exact matrix product states in the following sections.

3 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.

3.1 Time-dependent matrix product ansatz

The vector $\bm{\psi}(t)$ defined in (29) 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 (34)) using a matrix product ansatz

[TABLE]

The matrices $V_{\tau},W_{\tau}$, $\tau=0,\pm$, are operators acting on the auxiliary space ${\mathbb{C}}^{2}\otimes{\cal F}$, while the boundary vectors ${\langle l|}\in{\mathbb{C}}^{2}\otimes{\cal F}^{*}$ and ${|r\rangle}\in{\mathbb{C}}^{2}\otimes{\cal F}$ single out an element of contracted matrices which yields correct probability amplitudes $\psi^{t}_{\tau_{-t+1},\dots,\tau_{t}}$. For the matrices ($\epsilon=\pm$), we have the following explicit matrix representation:

[TABLE]

and for the boundary vectors we have

[TABLE]

where we have used the canonical basis vectors $\mathsf{e}_{j}$, $j=1,2$ for ${\mathbb{C}}^{2}$ and the $B$-eigenvector $\mathsf{v}=\mathsf{e}_{1}+\mathsf{e}_{2}$.

We now show that the matrix product expression (41) indeed provides the correct time-dependent 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 $\bm{\psi}(t)$ can be expressed as follows

[TABLE]

where $\bm{V}$ and $\bm{W}$ are the 3-component vectors encompassing the matrices

[TABLE]

The subscripts in (44) 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 (44) are equivalent to (41). The first set of algebraic relations concerns the commutation relations between the matrices and can be concisely written using the tensor space formalism

[TABLE]

that is to say in components

[TABLE]

The other two sets of algebraic relations concern the boundary vectors on the left and on the right and read

[TABLE]

or equivalently

[TABLE]

These relations can be checked by direct computation.

The algebraic relations ensure that the propagation equation is fulfilled

[TABLE]

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.

3.2 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 normalization $Z_{t}$ is defined as the sum of entries of the vector $\bm{\psi}(t)$. It can be checked from the explicit expression of the matrices (42) and boundary vectors (43) that $\bm{\psi}(t)$ is correctly normalized, i.e.

[TABLE]

where $T=(V_{0}+V_{+}+V_{-})(W_{0}+W_{+}+W_{-})$. In order to lighten the notations, we have introduced the parameters $\alpha=\alpha_{+}+\alpha_{-}$ and $\beta=\beta_{+}+\beta_{-}$, 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 $\epsilon$ at site $j=2k-t$, with $1\leq k\leq t$, can be easily expressed within the matrix product framework as

[TABLE]

This formula can be exactly evaluated using the explicit expression of the matrices and boundary vectors, see appendix A.2

[TABLE]

Note that the total particle density is constant

[TABLE]

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-\alpha$ for $t-j$ even.

The probability to observe a particle of species $\epsilon$ at site $j=2k+1-t$, with $0\leq k\leq t-1$, can be exactly computed using the matrix product expression, see appendix A.2

[TABLE]

Again, the total particle density has a very simple expression

[TABLE]

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\leq-t$ or $j>t$, it is straightforward to compute

[TABLE]

3.3 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

[TABLE]

for $j=2k-t$, i.e. in the case where $t-j$ is even. A similar expression also exists for $j=2k+1-t$

[TABLE]

The asymptotic behavior of the contour integral on the Euler scale $j=xt$, with $-1\leq x\leq 1$, and $t$ goes to infinity, can be studied using the saddle point analysis, which yields

[TABLE]

[TABLE]

At the junction of the two density profiles $x_{0}=\frac{\alpha-\beta}{\alpha+\beta}$, where we observe a discontinuity on the Euler scale, the profile exhibits a smooth transition on the diffusive scale $\frac{x-x_{0}}{t^{2}}$ [19, 20]. More precisely, for $j=x_{0}t+y\sqrt{t}$, we have

[TABLE]

where we recall the definition of the error function

[TABLE]

Similarly, we have

[TABLE]

The equations (62) and (64) are proven starting from the contour integral expressions (58) and (59), using the expansion $(\beta-\alpha z)^{-1}=\beta^{-1}\sum_{l=0}^{\infty}(z\alpha/\beta)^{l}$. A saddle point analysis reveals that the terms of the sum with non-vanishing contribution in the limit $t\to\infty$ are those for which $l$ is of the order $\sqrt{t}$, and yields the desired formulas.

It is instructive to characterize Fig. 3 from the hydrodynamic point of view, and interpret it in terms of the Riemann problem for a fluid with two locally conserved species (+ and - particles). As an outcome of temporal evolution of initially sharp profile $(\alpha_{+},\alpha_{-}|\beta_{+},\beta_{-})$ we have: two contact discontinuities of the type Eq.(20), Eq.(19), propagating along the light cone, and a domain wall (discontinuous on the Euler scale) between the two separable current carrying states of the type described in subsec. 2.2. In the following we shall denote a generic separable state from (14) with $\beta_{\pm}=\mu\alpha_{\pm}$ as $(\alpha_{\pm};\beta_{\pm})\equiv(\alpha_{\pm};\mu\alpha_{\pm})$, and an interface between two separable states $(\alpha_{\pm};\beta_{\pm})$ and $(\gamma_{\pm};\delta_{\pm})$ as $(\alpha_{\pm};\beta_{\pm}|\gamma_{\pm};\delta_{\pm})$.

With the above notations, on the left boundary of the light cone we have a contact discontinuity between a separable state with no current and a current-carrying separable state $(\alpha_{\pm};\alpha_{\pm}|\frac{\alpha}{\beta}\alpha_{\pm};\alpha_{\pm})$ of the type (20), propagating to the left with light velocity $-1$, that can be read off from Eqs.(60), (61). Likewise, at the right boundary of the light cone we have a contact discontinuity $(\beta_{\pm};\frac{\beta}{\alpha}\beta_{\pm}|\beta_{\pm};\beta_{\pm})$, propagating to the right with the velocity of light $+1$. As discussed in subsec. 2.2, these contact discontinuities correspond to interaction-free normal modes due to free propagation of the vacancies. Note that existence of two modes for vacancy propagation is due to even-odd staggered nature of the propagator.

Finally, there is a single diffusive domain wall (discontinuous on Euler scale), of the type $(\frac{\alpha}{\beta}\alpha_{\pm};\alpha_{\pm}|\beta_{\pm};\frac{\beta}{\alpha}\beta_{\pm})$. Its velocity can be computed from the mass conservation (22), (25) and it gives

[TABLE]

in accordance with exact analytic results (60), (61). In addition we find the characteristic velocities of the diffusive modes at both sides of the shock to be $c_{1}^{L}=c_{1}^{R}=\frac{\alpha-\beta}{\alpha+\beta}=v_{sh}$. According to the Lax criterium, the shock on Fig. 3 is marginally stable (a diffusive discontinuity is parallel to the mode velocities on both sides of the discontinuity).

4 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.

4.1 Time-dependent matrix product ansatz

The vector $\bm{\psi}(t)$ defined in (32) 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 (34)) using a matrix product ansatz

[TABLE]

Note that the ansatz is a bit more complicated than previously because it involves a boundary vector ${\langle l|}_{\tau}$ which depends on the content of the site $-t+1$. This dependence is reminiscent of the defect in the initial condition and could be equivalently encoded using a defect matrix (i.e different from $V_{\tau},W_{\tau}$) on the site $-t+1$. The matrices $V_{\tau},W_{\tau}$, $\tau=0,\pm$, are again operators acting in the auxiliary space ${\mathbb{C}}^{2}\otimes{\cal F}$. The boundary vectors ${\langle l|}_{\tau}\in{\mathbb{C}}^{2}\otimes{\cal F}^{*}$ and ${|r\rangle}\in{\mathbb{C}}^{2}\otimes{\cal 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

[TABLE]

Note that $\mathsf{P}$ is a projector ($\mathsf{P}^{2}=\mathsf{P}$) obeying $V_{\tau}\,\mathsf{P}=V_{\tau}$ and $W_{\tau}\,\mathsf{P}=W_{\tau}$, $\tau=\pm\,,0$. For the boundary vectors, we choose the following form

[TABLE]

where we have used a coherent state in ${\langle l|}_{0}$, and canonical states in ${\langle l|}_{\pm}$ and ${|r\rangle}$. Note that ${\langle l|}_{\tau}\,\mathsf{P}={\langle l|}_{\tau}$, $\tau=\pm\,,0$.

In order to efficiently show that the matrix product expression (66) provides the correct time-dependent probability distribution, we reformulate the matrix product expression using the tensor space formalism

[TABLE]

where $\bm{V}$, $\bm{W}$ and $\bm{{\langle l|}}$ are the 3-component vectors

[TABLE]

The subscripts in (69) 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 (46). On the left boundary, the following condition should be satisfied

[TABLE]

that is to say

[TABLE]

On the right boundary, we have

[TABLE]

and finally, the boundary vectors should also satisfy the ‘initial condition’

[TABLE]

The algebraic relations ensure that the propagation equation is fulfilled

[TABLE]

Once again the matrix product structure will prove to be very efficient to analytically compute physical quantities, such as the density profiles.

4.2 Exact expression of physical observables

The first step toward the exact computation of the time-dependent density profiles is to check that the vector $\bm{\psi}(t)$, defining the probability distribution in the propagation cone, is correctly normalized. The normalization $Z_{t}$ of the vector $\bm{\psi}(t)$ is defined as the sum of its components and can be written using the matrix product formalism as follows

[TABLE]

where $T=(W_{0}+W_{+}+W_{-})(V_{0}+V_{+}+V_{-})$. It can be evaluated, see appendix B.1, using the explicit form of the matrices and of the boundary vectors

[TABLE]

where we have introduced

[TABLE]

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 $\epsilon$ at site $j=2k-t$, with $1\leq k\leq t$, can be easily expressed within the matrix product framework as

[TABLE]

This formula can be exactly evaluated using the explicit expression of the matrices and boundary vectors, see appendix B.2

[TABLE]

The total particle density takes the very simple expression

[TABLE]

The probability to observe a particle of species $\epsilon$ at site $j=2k+1-t$, with $1\leq k\leq t-1$, can be exactly computed using the matrix product expression, see appendix B.2

[TABLE]

Once again, the total particle density is simply

[TABLE]

The density on the site $j=-t+1$ (i.e. for $k=0$) takes a specific expression, see appendix B.2

[TABLE]

In this case the total particle density reads

[TABLE]

4.3 Hydrodynamic limit

An equivalent alternative expression for the density is given by

[TABLE]

Using the saddle point analysis we can again study the asymptotic behavior of the contour integral in the regime where $j=xt$, with $-1\leq x\leq 1$, and $t$ goes to infinity. The contribution of the contour integral is exponentially small when $x\neq 0$, and one gets

[TABLE]

Note that at the site $j=-t+1$, we can obtain an explicit expression, which in the long time limit reads

[TABLE]

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=\frac{y}{\sqrt{t}}$, we obtain

[TABLE]

The asymptotic behavior of n_{\epsilon}^{\mathsf{o}}\big{(}y\sqrt{t},t\big{)} 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).

5 Stochastic boundary driving

5.1 Matrix product expression of the stationary state

We recall that we are using the short-hand notations

[TABLE]

The main result of the present subsection is an exact expression of the components of the stationary states $\bm{p}$, defined by the equations ${\mathbb{U}}\bm{p}=\bm{p}$ (where the operator ${\mathbb{U}}$ has been defined in the subsection 2.3). We have

[TABLE]

where the matrices $V_{\tau},W_{\tau}\in\mbox{End}({\mathbb{C}}^{2}\otimes{\cal F}\otimes{\cal 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 ${\cal F}$. The matrix product ansatz is constructed by employing matrices $A$ and $B$, which were introduced in the preceding section (40):

[TABLE]

The vectors ${\langle l|}_{\tau}\in{\mathbb{C}}^{2}\otimes{\cal F}^{*}\otimes{\cal F}^{*}$ read explicitly

[TABLE]

The vectors ${|rr\rangle}_{\tau,\tau^{\prime}}\in{\mathbb{C}}^{2}\otimes{\cal F}\otimes{\cal F}$ are expressed using the $B$-eigenvector $\mathsf{v}=\mathsf{e}_{1}+\mathsf{e}_{2}$

[TABLE]

In order to prove the matrix product expression of the stationary state it will be useful to reformulate it using the tensor space formalism

[TABLE]

where

[TABLE]

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 $\bm{p}$ and $\bm{p^{\prime}}$ satisfying the following relations

[TABLE]

These relations imply the following eigenvector relation

[TABLE]

We already provided the expression of the stationary state $\bm{p}$. We now give the expression of the vector $\bm{p^{\prime}}$. For this purpose we need to introduce the boundary vectors

[TABLE]

where

[TABLE]

and

[TABLE]

Building on those, we define

[TABLE]

The relations (98) are proven using local exchange relations in the bulk

[TABLE]

on the right boundary

[TABLE]

and on the left boundary

[TABLE]

By direct computation it is possible to check that all of the above relations are satisfied by the ansatz (93), (94) and (95).

5.2 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

[TABLE]

where $T=(V_{0}+V_{+}+V_{-})(W_{0}+W_{+}+W_{-})$. We show in appendix C.1 that it is given by the following exact formula

[TABLE]

The mean particle current associated to the species $\epsilon=\pm$ is defined by the mean number of particles of species $\epsilon$ 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

[TABLE]

The above expression can be evaluated exactly (see appendix C.2)

[TABLE]

From this result we can deduce that the total particle current takes a very simple form

[TABLE]

When $\alpha=\beta$, the current is given by the appropriate limit of the general expression

[TABLE]

The mean particle density associated to the species $\epsilon=\pm$ at site $i$ is defined by the probability to observe a particle $\epsilon$ 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 $\bm{p}$ and $\bm{p^{\prime}}$ respectively. The density at site $i=2k+1$ is thus given in the matrix product formalism by

[TABLE]

where $\tilde{Z}_{L}$ is the normalization (i.e. the sum of the components) of $\bm{p^{\prime}}$. Once again the matrix product expression can be computed explicitly (see appendix C.3)

[TABLE]

We deduce that the total particle density takes a very simple form

[TABLE]

Note that we have an equivalent expression of the particle density (114)

[TABLE]

which allows us to easily take the limit $\beta\to\alpha$, and obtain an expression of the particle density in the case $\alpha=\beta$

[TABLE]

Note that in this case the density profile is a linear profile interpolating between the densities of the left and right reservoirs $\alpha_{\epsilon}$ and $\beta_{\epsilon}$.

5.3 Hydrodynamic limit and phase transitions

In this subsection we investigate the properties of the physical quantities in the large system size limit $L\to\infty$. 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 $\alpha<\beta$.

The particle current has a finite large system size limit

[TABLE]

The limit of the particle density is

[TABLE]

It has a constant value $\frac{1}{2}\frac{\beta_{\epsilon}}{\beta}(\alpha+\beta)$, fixed by the right reservoir in the bulk, and an exponential tail near the left boundary.

The left reservoir phase, for $\beta<\alpha$.

The particle current has a finite large system size limit

[TABLE]

The limit of the particle density is

[TABLE]

It has a constant value $\frac{1}{2}\frac{\alpha_{\epsilon}}{\alpha}(\alpha+\beta)$, fixed by the left reservoir, and an exponential tail near the right boundary.

The diffusive phase, for $\alpha=\beta$.

In this phase the particle current decays as $1/L$ in the large system size limit. More precisely we have

[TABLE]

The limit of the particle density is

[TABLE]

In this phase the Fick law (which is a signature of the diffusive behavior) is satisfied

[TABLE]

where ${\cal D}=\frac{1-\alpha}{\alpha}$ is the diffusion constant, in accordance with the value computed in [19, 20].

6 Discussion and 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.

Although in this paper we focused on a model with deterministic dynamics in the bulk, we believe that very similar methods can be employed also to solve certain interesting stochastic deformations of the model. For example, in Ref. [19], a stochastic variant of our lattice dynamics has been proposed where the local propagator (6) is replaced by a Markov matrix

[TABLE]

where $\Gamma\in[0,1]$ is a probability of exchange of $+$ and $-$ particles upon their collision (scattering). Remarkably, even though $U_{\Gamma}$ in general no longer satisfies the braid relation (7), we find that certain nontrivial dynamical [19] and non-equilibrium steady state problems can still be exactly solved. In particular, we have empirical evidence that the stochastic deformation (125) of the steady state driven by stochastic boundaries can be written by a matrix product ansatz of exactly the same complexity (characterised by bond dimensions, or Schmidt ranks of all bi-partitions) as in the deterministic case ($\Gamma=0$) studied in sect. 5.1. However, its precise analytic structure shall be left for future work.

Acknowledgements

We warmly thank Katja Klobas for useful discussions.

Funding information

Support from the Advanced grant 694544-OMNES of the European Research Council (ERC), and the Program P1-0402 of Slovenian Research Agency (ARRS) are gratefully acknowledged. VP acknowledges support from Deutsche Forschungsgemeinschaft, and from interdisciplinary UoC Forum “Classical and quantum dynamics of interacting particle systems”.

Appendix A Details of the computations for the inhomogeneous quench

A.1 Computation of the normalisation

The normalisation of the time-dependent matrix product state is defined as

[TABLE]

To evaluate explicitly this quantity, we first use the fact that

[TABLE]

which can be readily checked using the explicit expression of the matrices (42). The normalisation can then be rewritten as

[TABLE]

The boundary vectors ${\langle l|}$ and ${|r\rangle}$ (see explicit expressions in (43)) fulfill the following eigenvalue relations

[TABLE]

Those relations immediately yield the expression

[TABLE]

A.2 Computation of the densities

The probability to observe a particle of species $\epsilon$ at site $j=2k-t$, with $1\leq k\leq t$, is given by the matrix product expression

[TABLE]

where $T=(W_{0}+W_{+}+W_{-})(V_{0}+V_{+}+V_{-})$. Using the commutation property (127) and the eigenvectors relations (129) we can simplify the expression

[TABLE]

The next step is to use the closure relation in the Fock space

[TABLE]

and insert it to the left of $W_{\epsilon}$. Note that the integration contour includes the pole at $z=0$ and excludes all others. Then using the fact that the coherent states ${|\underline{z}\rangle}$ and ${\langle\underline{1/z}|}$ are respectively eigenstates of the annihilation and creation operators $\mathfrak{a}$ and $\mathfrak{a}^{\dagger}$, we obtain

[TABLE]

The integrand can be evaluated using the triangular structure of matrices $A$ and $B$

[TABLE]

The residue at $z=0$ corresponding to the first of the two previous terms is equal to

[TABLE]

whereas the residue at $z=0$ corresponding to the second term is equal to

[TABLE]

Gathering these results together we end up with the following expression

[TABLE]

Finally, using the fact that $\alpha_{0}+\alpha=1$ and $\beta_{0}+\beta=1$, the expression reduces to

[TABLE]

Note that we also have an alternative expression for the density, if we do not evaluate explicitly the residue (136)

[TABLE]

where the contour of the integral encircles the pole at $z=0$ and excludes the pole at $z=\beta/\alpha$. The probability to observe a particle of species $\epsilon$ at site $j=2k+1-t$, with $0\leq k\leq t-1$, is given by the matrix product expression

[TABLE]

A very similar computation yields

[TABLE]

Using the fact that $\alpha_{0}+\alpha=1$ and $\beta_{0}+\beta=1$, the expression reduces to

[TABLE]

As before, we also have a contour integral expression

[TABLE]

Appendix 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

[TABLE]

where $T=(W_{0}+W_{+}+W_{-})(V_{0}+V_{+}+V_{-})$. To evaluate this quantity explicitly we first use the fact that

[TABLE]

which can be readily checked using the explicit expression of the matrices (67). The normalization can then be rewritten as

[TABLE]

The boundary vectors ${\langle l|}_{\tau}$ and ${|r\rangle}$ (see explicit expressions in (68)) fulfill the following eigenvalue relations

[TABLE]

We have also the scalar product expression

[TABLE]

Those relations lead immediately to the expression

[TABLE]

B.2 Computation of the densities

The probability to observe a particle of species $\epsilon$ at site $j=2k-t$, with $1\leq k\leq t$, is expressed within the matrix product framework as

[TABLE]

Using the commutation property (146) and the eigenvectors properties (148) the expression is simplified into

[TABLE]

Then using the explicit expression of the matrices (67) and of the boundary vectors (68) and the properties $V_{\tau}\,\mathsf{P}=V_{\tau}$, $W_{\tau}\,\mathsf{P}=W_{\tau}$ and ${\langle l|}_{\tau}\,\mathsf{P}={\langle l|}_{\tau}$, $\tau=\pm\,,0$, we can get rid of the projector $\mathsf{P}$ in the matrix product and we obtain

[TABLE]

The next step is to use the closure relation in the Fock space

[TABLE]

and insert it once between222you can indeed observe in equation (153) that the operators $\mathfrak{a}$ are all located on the left side and the operators $\mathfrak{a}^{\dagger}$ are all located on the right side of the expression the operators $\mathfrak{a}$ and $\mathfrak{a}^{\dagger}$. Note that the integration contour includes the pole at $z=0$ and excludes all others. Then using the fact that the coherent states ${|\underline{z}\rangle}$ and ${\langle\underline{1/z}|}$ are respectively eigenstates of the annihilation and creation operators $\mathfrak{a}$ and $\mathfrak{a}^{\dagger}$, we obtain

[TABLE]

The integrand can be further simplified using the triangular structure of the matrices. It reads

[TABLE]

A straightforward computation yields the following results. The residue at $z=0$ corresponding to the third term is equal to [math]. The residue corresponding to the second term at $z=0$ is equal to $\rho_{\epsilon}(\lambda_{0}+\lambda)(\rho_{0}+\rho)^{t-1}$, whereas the one corresponding to the first term is given by

[TABLE]

Gathering the results together we obtain an expression of the density

[TABLE]

Using the fact that $\rho_{0}+\rho=1$ and $\lambda_{0}+\lambda=1$ the expression reduces to

[TABLE]

Note that we have also an alternative expression of the density if we do not evaluate explicitly the residue associated to the first term in (156)

[TABLE]

where the contour of the integral encircles the pole at $z=0$. The probability to observe a particle of species $\epsilon$ at site $j=2k+1-t$, with $1\leq k\leq t-1$, is given by the matrix product expression

[TABLE]

A very similar computation yields to

[TABLE]

which reduces to

[TABLE]

If we do not evaluate explicitly the residue, we obtain an alternative expression of the densities,

[TABLE]

Finally we have the particular case of the density of species $\epsilon$ at site $j=-t+1$. The corresponding matrix product expression is

[TABLE]

Using the commutation property (146), the right eigenvector property in (148) we obtain

[TABLE]

Then using recursively the relations

[TABLE]

we get

[TABLE]

The scalar product relations

[TABLE]

yield the final expression

[TABLE]

Taking into account the fact that $\rho_{0}+\rho=1$ and $\lambda_{0}+\lambda=1$, the formula reduces to

[TABLE]

Appendix C Details of the computations for the boundary driven model

C.1 Computation of the normalisation

The normalisation of the stationary distribution is defined as

[TABLE]

The first step to compute this quantity is to remark the commutation property

[TABLE]

which can be directly checked using the explicit representation (93). We can thus rewrite

[TABLE]

We have the following expressions

[TABLE]

Using these expressions it is straightforward to see that the actions of the coherent states ${\langle\underline{\beta/\beta_{0}}|}$ (in the third tensor component) and ${|\underline{\alpha/\alpha_{0}}\rangle}$ (in the second tensor component) can be evaluated because of the relations (38). We thus have

[TABLE]

and

[TABLE]

Using the explicit expression

[TABLE]

we can compute

[TABLE]

and

[TABLE]

Gathering all these results together, the normalisation can be explicitly evaluated

[TABLE]

Taking into account that $\alpha+\alpha_{0}=1$ and $\beta+\beta_{0}=1$, the expression reduces to

[TABLE]

C.2 Computation of the particle currents

The particle current associated to the species $\epsilon=\pm$ is expressed as (using again the commutation property (173))

[TABLE]

We have the relation

[TABLE]

Similarly to what has been done for the computation of the normalisation, we can evaluate the action of the coherent states ${|\underline{\alpha/\alpha_{0}}\rangle}$ and ${\langle\underline{\beta/\beta_{0}}|}$. It allows one to realize that the contributions of $\alpha_{\epsilon}\alpha(\beta+\beta_{0})\mathsf{e}_{1}^{t}\otimes{\langle 0|}\otimes{\langle\underline{\beta/\beta_{0}}|}$ and of $-\alpha_{0}\alpha_{\epsilon}(\beta+\beta_{0})\mathsf{e}_{1}^{t}\otimes{\langle 1|}\otimes{\langle\underline{\beta/\beta_{0}}|}$ (appearing in (184)) in the computation of the current cancel with each other because ${\langle 0|}\cdot{|\underline{\alpha/\alpha_{0}}\rangle}=1$ and ${\langle 1|}\cdot{|\underline{\alpha/\alpha_{0}}\rangle}=\alpha/\alpha_{0}$. Concerning the remaining terms, we have

[TABLE]

and

[TABLE]

The contribution of the term $\alpha_{\epsilon}\beta_{0}(\alpha+\alpha_{0})\mathsf{e}_{1}^{t}\otimes{\langle 1|}\otimes{\langle 0|}$ of (184) is a bit more involved to compute because we do not have the coherent state ${\langle\underline{\beta/\beta_{0}}|}$ anymore. We are going to show that

[TABLE]

Using (184), (185), (186) and (187) it is then straightforward to show that

[TABLE]

Taking into account that $\alpha+\alpha_{0}=1$ and $\beta+\beta_{0}=1$, the expression reduces to

[TABLE]

We now come back to the proof of equation (187). The first step is to observe that the action of the coherent state ${|\underline{\alpha/\alpha_{0}}\rangle}$ in the second tensor component of the right boundary vectors can be evaluated, similarly as in the previous computations. The action of the vector ${\langle 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 (175)), so that $\mathsf{e}_{1}^{t}$ is a left eigenvector of the remaining terms ($\mathsf{e}_{1}^{t}$ is a left eigenvector of the matrix $A$). Applying all these manipulations we end up with

[TABLE]

The last point left to complete the proof is to show the following relation

[TABLE]

Using the triangular structure of the matrices $A$ and $B$ we can evaluate the matrix product in the finite dimensional auxiliary space as follows

[TABLE]

To compute the expression above, we introduce the ket and bra vectors

[TABLE]

where $\gamma=\sqrt{\frac{\alpha_{0}\beta_{0}}{\alpha\beta}}$ and $\delta=\sqrt{\frac{\alpha\beta_{0}}{\alpha_{0}\beta}}$. Note that they can be conveniently expressed using coherent states

[TABLE]

and similarly for ${\langle w(z)|}$. They are the right and left eigenvectors of the following operator

[TABLE]

They satisfy also the following closure relation

[TABLE]

where the integration contour includes poles at $z=0$ and at $z=\delta$ and excludes any other. We insert the closure relation in (192) and then use the eigenvalue relations (195) and the fact that

[TABLE]

The summation in (192) is then simply a geometric summation in the integrand. After simplifications we are left with

[TABLE]

The residue at $z=\delta$ of the integrand can be easily computed (it is a simple pole) and is equal to

[TABLE]

The computation of the residue at $z=0$ is a bit more involved due to multiple order poles coming from the expansion of $\frac{1}{z}\left(1+\frac{1}{\gamma z}\right)^{n}$. The residue at $z=0$ is

[TABLE]

The $l+1$ order derivative can be computed using the simple elements decomposition

[TABLE]

and the fact that

[TABLE]

Gathering all terms together yields the following expression of the residue at $z=0$

[TABLE]

which implies that the contour integral in (198) is equal to $\left(1+\frac{\alpha}{\alpha_{0}}\right)^{n}$ and concludes the proof of the formula (191).

C.3 Computation of the densities

The mean particle density associated to the species $\epsilon=\pm$ at site $i$ is defined by the probability to observe a particle of species $\epsilon$ at site $i$ in the stationary state. We need to average over the two steps (odd and even) of the dynamics: in the stationary state the system is half of the time described by the probability vector $\bm{p}$ and the other half by the probability vector $\bm{p^{\prime}}$. The density at site $i=2k+1$ is thus given by

[TABLE]

where $T=(W_{0}+W_{+}+W_{-})(V_{0}+V_{+}+V_{-})$ and $\tilde{Z}_{L}$ is the normalisation of $\bm{p^{\prime}}$. Using the explicit expression of the matrices (93) and boundary vectors (94), (95), (102) and (101) we can use very similar techniques that have been used to compute the particle currents (partial action of the coherent states and the closure relation) to show that

[TABLE]

and

[TABLE]

Summing both contribution we end up with

[TABLE]

Taking into account that $\alpha+\alpha_{0}=1$ and $\beta+\beta_{0}=1$, the expression reduces to

[TABLE]