Exact matrix product decay modes of a boundary driven cellular automaton

TL;DR

This paper investigates the integrability and spectral properties of a boundary-driven cellular automaton, providing explicit decay modes, Bethe-like equations, and insights into the spectrum's structure and thermodynamic behavior.

Contribution

It introduces an exact matrix product ansatz for decay modes and conjectures Bethe-like equations for all spectral orbitals, advancing understanding of boundary-driven automaton dynamics.

Findings

Explicit decay modes constructed for the automaton

Spectrum separates into orbitals with Bethe-like equations

Thermodynamic properties and spectral gap scaling analyzed

Abstract

We study integrability properties of a reversible deterministic cellular automaton (the rule 54 of [Bobenko et al., Commun. Math. Phys. 158, 127 (1993)]) and present a bulk algebraic relation and its inhomogeneous extension which allow for an explicit construction of Liouvillian decay modes for two distinct families of stochastic boundary driving. The spectrum of the many-body stochastic matrix defining the time propagation is found to separate into sets, which we call orbitals, and the eigenvalues in each orbital are found to obey a distinct set of Bethe-like equations. We construct the decay modes in the first orbital (containing the leading decay mode) in terms of an exact inhomogeneous matrix product ansatz, study the thermodynamic properties of the spectrum and the scaling of its gap, and provide a conjecture for the Bethe-like equations for all the orbitals and their degeneracy.

Equations

$$s_{\mathrm{S}}=\chi (s_{\mathrm{W}},s_{\mathrm{N}},s_{\mathrm{E}})=s_{\mathrm{N}}+s_{\mathrm{W}}+s_{\mathrm{E}}+s_{\mathrm{W}}{s}_{\mathrm{E}}(\mathrm{mod}2),s_{\mathrm{S}},s_{\mathrm{N}},s_{\mathrm{W}},s_{\mathrm{E}}\in {\mathbb{Z}}_2.$$

$$P_{(s,s^{\prime },s^{\prime \prime }),(t,t^{\prime },t^{\prime \prime })}={\delta }_{s,t}{\delta }_{s^{\prime },\chi (t,t^{\prime },t^{\prime \prime })}{\delta }_{s^{\prime \prime },t^{\prime \prime }},$$

$$P=\left(\begin{array}{cccccccc}1& & & & & & & \\ & & & 1& & & & \\ & & 1& & & & & \\ & 1& & & & & & \\ & & & & & & 1& \\ & & & & & & & 1\\ & & & & 1& & & \\ & & & & & 1& & \end{array}\right),$$

$$\mathbf{p}(t)=U^t\mathbf{p}(0),$$

$$U=U_{\mathrm{o}}{U}_{\mathrm{e}},$$

$$U_{\mathrm{e}}=P_{123}{P}_{345}\cdots P_{n-3,n-2,n-1}{P}_{n-1,n}^{\mathrm{R}},$$

$$U_{\mathrm{o}}=P_{12}^{\mathrm{L}}{P}_{234}{P}_{456}\cdots P_{n-2,n-1,n}.$$

$$P_{12}^{\mathrm{L}}=P^{\mathrm{L}}\otimes \text{mathbbm}{1}_{2^{n-2}},P_{n-1,n}^{\mathrm{R}}=\text{mathbbm}{1}_{2^{n-2}}\otimes P^{\mathrm{R}}.$$

$$[P_{k-1,k,k+1},P_{k^{\prime }-1,k^{\prime },k^{\prime }+1}]=0,\mathrm{if}|k-k^{\prime }|\ge 2.$$

$$[P_{12}^{\mathrm{L}},P_{234}]=0,[P_{123},P_{34}^{\mathrm{R}}]=0,$$

$$U\mathbf{p}=\Lambda \mathbf{p},$$

$$U_{\mathrm{e}}\mathbf{p}={\Lambda }_{\mathrm{L}}{\mathbf{p}}^{\prime },U_{\mathrm{o}}{\mathbf{p}}^{\prime }={\Lambda }_{\mathrm{R}}\mathbf{p},\Lambda ={\Lambda }_{\mathrm{R}}{\Lambda }_{\mathrm{L}}.$$

$$\mathbf{p}(t)={\mathbf{p}}_0+\sum _{j\ge 1}{c}_j{\Lambda }_j^t{\mathbf{p}}_j,$$

$$\displaystyle P^{\rm{L}}=\left(\begin{array}[]{cccc}\alpha&0&\alpha&0\\ 0&\beta&0&\beta\\ 1-\alpha&0&1-\alpha&0\\ 0&1-\beta&0&1-\beta\\ \end{array}\right),\quad P^{\rm{R}}=\left(\begin{array}[]{cccc}\gamma&\gamma&0&0\\ 1-\gamma&1-\gamma&0&0\\ 0&0&\delta&\delta\\ 0&0&1-\delta&1-\delta\\ \end{array}\right),$$

$$\displaystyle P^{\rm{L}}=\left(\begin{array}[]{cccc}\frac{1}{2}&0&\frac{1}{2}&0\\ 0&1-\alpha&0&\beta\\ \frac{1}{2}&0&\frac{1}{2}&0\\ 0&\alpha&0&1-\beta\\ \end{array}\right),\qquad P^{\rm{R}}=\left(\begin{array}[]{cccc}\frac{1}{2}&\frac{1}{2}&0&0\\ \frac{1}{2}&\frac{1}{2}&0&0\\ 0&0&1-\gamma&\delta\\ 0&0&\gamma&1-\delta\\ \end{array}\right),$$

$$\mu :=\alpha -\beta ,\nu :=\beta -1,\sigma :=\gamma -\delta ,\rho :=\delta -1.$$

$$\mathbf{W}=\left(\begin{array}{c}{W}_0(\xi ,\omega )\\ W_1(\xi ,\omega )\end{array}\right),$$

$$W_{0}(\xi,\omega)=\left(\begin{array}[]{cccc}1&1&0&0\\ 0&0&0&0\\ \xi&\xi&0&0\\ 0&0&0&0\\ \end{array}\right),\qquad W_{1}(\xi,\omega)=\left(\begin{array}[]{cccc}0&0&0&0\\ 0&0&\xi&1\\ 0&0&0&0\\ 0&0&1&\omega\\ \end{array}\right),$$

$${\mathbf{W}}^{\prime }(\xi ,\omega )=\mathbf{W}(\omega ,\xi ).$$

$$P_{123}{\mathbf{W}}_{1}S{\mathbf{W}}_2{\mathbf{W}}_3^{\prime }={\mathbf{W}}_1{\mathbf{W}}_2^{\prime }{\mathbf{W}}_{3}S,$$

$$S=\mathbbm{1}\otimes\sigma^{\rm{x}}=\left(\begin{array}[]{cccc}0&1&0&0\\ 1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ \end{array}\right),$$

$$P_{123}{\mathbf{W}}_1^{\prime }{\mathbf{W}}_2{\mathbf{W}}_3^{\prime }S={\mathbf{W}}_1^{\prime }S{\mathbf{W}}_2^{\prime }{\mathbf{W}}_3.$$

$$\mathbf{p}=⟨{\mathbf{l}}_1|{\mathbf{W}}_2{\mathbf{W}}_3^{\prime }{\mathbf{W}}_4{\mathbf{W}}_5^{\prime }\cdots {\mathbf{W}}_{n-3}^{\prime }{\mathbf{W}}_{n-2}|{\mathbf{r}}_{n-1,n}⟩,$$

$${\mathbf{p}}^{\prime }=⟨{\mathbf{l}}_{12}^{\prime }|{\mathbf{W}}_3{\mathbf{W}}_4^{\prime }{\mathbf{W}}_5{\mathbf{W}}_6^{\prime }\cdots {\mathbf{W}}_{n-2}^{\prime }{\mathbf{W}}_{n-1}|{\mathbf{r}}_n^{\prime }⟩.$$

$$P_{12}^{\mathrm{L}}⟨{\mathbf{l}}_{12}^{\prime }|={\lambda }_{\mathrm{L}}⟨{\mathbf{l}}_1|{\mathbf{W}}_{2}S,$$

$$P_{123}⟨{\mathbf{l}}_1|{\mathbf{W}}_2{\mathbf{W}}_3^{\prime }=⟨{\mathbf{l}}_{12}^{\prime }|{\mathbf{W}}_{3}S,$$

$$P_{12}^R|{\mathbf{r}}_{12}⟩={\mathbf{W}}_1^{\prime }S|{\mathbf{r}}_2^{\prime }⟩,$$

$$P_{123}{\mathbf{W}}_1^{\prime }{\mathbf{W}}_2|{\mathbf{r}}_3^{\prime }⟩={\lambda }_{\mathrm{R}}{\mathbf{W}}_1^{\prime }S|{\mathbf{r}}_{23}⟩.$$

$$\xi =\frac{(\alpha +\beta -1)-\beta {\lambda }_{\mathrm{L}}}{(\beta -1){\lambda }_L^2},$$

$$\omega =\frac{{\lambda }_{\mathrm{L}}(\alpha -{\lambda }_{\mathrm{L}})}{(\beta -1)},$$

$$⟨l_0|=\left(0,\frac{1}{{\lambda }_{\mathrm{L}}},\frac{\alpha (\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{\beta {\lambda }_{\mathrm{L}}^2(\alpha -{\lambda }_{\mathrm{L}})},1\right),$$

$$⟨l_1|=\left(0,\frac{1-\beta }{\beta {\lambda }_{\mathrm{L}}},\frac{(1-\alpha )(\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{\beta {\lambda }_{\mathrm{L}}^2(\alpha -{\lambda }_{\mathrm{L}})},\frac{1}{\beta }-1\right),$$

$$⟨l_{0,0}^{\prime }|=\left(\frac{\alpha (\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{(\beta -1)\beta {\lambda }_{\mathrm{L}}},\frac{\alpha +\beta -1-{\lambda }_{\mathrm{L}}}{\beta -1},0,0\right),$$

$$⟨l_{0,1}^{\prime }|=\left(0,0,\frac{{\lambda }_{\mathrm{L}}(\alpha -{\lambda }_{\mathrm{L}})(\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{(\beta -1)(\alpha +\beta -1-\beta {\lambda }_{\mathrm{L}})},\frac{\alpha (\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{(\beta -1)\beta {\lambda }_{\mathrm{L}}}\right),$$

$$⟨l_{1,0}^{\prime }|=\left(\frac{(\alpha -{\lambda }_{\mathrm{L}})(\alpha +\beta -1-{\lambda }_{\mathrm{L}})}{(1-\beta )\beta {\lambda }_{\mathrm{L}}},-\frac{\alpha +\beta -1-{\lambda }_{\mathrm{L}}}{\beta },0,0\right),$$

Full text

1 Introduction

Understanding the emergence of laws governing macroscopic physical phenomena, such as transport and relaxation, from deterministic and reversible microscopic dynamics is one of the most prominent fundamental problems of statistical mechanics. In this context, an important setup consists of driving a finite many-body system, say a one dimensional lattice with local interactions, with a pair of macroscopic (infinite) reservoirs attached, or coupled to the system’s ends (boundaries). Infinite reservoirs can typically be replaced by stochastic forces acting on boundary degrees of freedom of the system, so we are speaking of boundary driven deterministic dynamics. Several non-trivial exactly solvable examples of current carrying non-equilibrium steady states (NESS) of this type of dynamics have been recently found in the realm integrable quantum lattice models [1], however all attempts of exact constructions of dynamical modes of relaxation have failed so far.

Recently, NESS of a boundary driven reversible cellular automaton, which can be understood as a simple caricature of deterministic interacting dynamics, has been found [2] and its construction exhibits certain interesting algebraic properties. The cellular automaton is the Rule 54 of Bobenko et al [3], which is a two-state fully deterministic, reversible many-body interacting system and admits non-trivially scattering solitons.

The rule is given by a deterministic local mapping on a diamond-shaped plaquette $\chi:\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\to\mathbb{Z}_{2}$. A south site $s_{\rm S}$ is determined by a north, west and east sites

[TABLE]

Time runs in the north to south direction (see Fig. 1) and defines a simple interacting dynamics over a $1+1$ dimensional lattice $s_{x,t+1}=\chi(s_{x-1,t},s_{x,t-1},s_{x+1,t})$, where only lattice sites $(x,t)$ of fixed parity of $x\pm t$ are considered.

We shall now define dynamics over a finite chain of even number of sites $n$ with the initial data given by a configuration along a saw $(s_{1},s_{2},\ldots,s_{n})\equiv(s_{1,t+1},s_{2,t},s_{3,t+1},s_{4,t},\ldots,s_{n-1,t+1},s_{n,t})$ (see Fig. 2), which can be given as a composition of even site updates $s_{2y,t+1}=\chi(s_{2y-1,t},s_{2y,t-1},s_{2y+1,t})$ and odd site updates $s_{2y+1,t+2}=\chi(s_{2y,t+1},s_{2y+1,t},s_{2y+2,t+1})$. The dynamics is fully deterministic, except for the sites near the boundaries, where we shall prescribe appropriate local Markov stochastic processes by which we drive the model out of equilibrium.

We thus proceed to formulate the time evolution of the full probability distribution $p_{(s_{1},s_{2},\ldots,s_{n})}$, which we call a state.111This notion of the state (as a macro-state) should be distinguished from a binary state of an automaton. Since the meaning of the term should be clear from the context we use the same word for the two concepts. The state space is a convex subset of a vector space ${\cal S}=\mathbb{R}^{\cal C}=(\mathbb{R}^{2})^{\otimes n}$ of probability distributions over configurations $\mathbf{p}=(p_{0},p_{1},\ldots,p_{2^{n}-1})\in{\cal S}$ satisfying the non-negativity and normalization conditions, $p_{s}\geq 0$, $\sum_{s=0}^{2^{n}-1}p_{s}=1$. Here, $s$ is a binary coded configuration $s=\sum_{k=1}^{n}2^{n-k}s_{k}$. The local rule 54 can then be given in terms of a three-site $2^{3}\times 2^{3}$ permutation matrix $P$

[TABLE]

or

[TABLE]

such that $P^{2}=\mathbbm{1}$. The local update rule is in turn embedded into ${\rm End}({\cal S})$ as $P_{k,k+1,k+2}=\mathbbm{1}_{2^{k-1}}\otimes P\otimes\mathbbm{1}_{2^{n-k-2}}$ acting on any triple of neighboring sites $k,k+1,k+2$. The time evolution of the state vector $\mathbf{p}(t)\in{\cal S}$ starting from some initial state $\mathbf{p}(0)$ is written as

[TABLE]

where

[TABLE]

is the one-step propagator that is factored in terms of two temporal layers which generate staggered dynamics for, respectively, even and odd sites

[TABLE]

The boundary propagators

[TABLE]

are given in terms of $4\times 4$ stochastic matrices222By definition, stochastic matrices have non-negative elements which in each column sum to 1. $P^{\rm L}$ and $P^{\rm R}$ (to be specified later), which in turn imply that the full $2^{n}\times 2^{n}$ propagator $U$ itself is a stochastic matrix and thus conserves total probability during the time evolution. This dynamics which is bulk-deterministic and boundary-stochastic should be contrasted with related, though distinct, discrete time asymmetric exclusion process models [4, 5, 6], which feature both stochastic bulk dynamics as well as stochastic driving.

We note that $P_{k-1,k,k+1}$ changes only the site $k$, conditioned on the states of the sites $k-1$ and $k+1$, so the local propagators commute if at least two sites apart

[TABLE]

Furthermore we shall require that also boundary stochastic matrices commute with the neighbouring bulk propagators

[TABLE]

so the order of factors in either $U_{\rm e}$ or $U_{\rm o}$, Eqs. (5,6), is irrelevant.

The main objective of this paper is an exact solution of an eigenvalue equation for the Markov propagator

[TABLE]

which can be conveniently split into a pair of equations for an eigenstate at even and odd time layers with the eigenvalue $\Lambda$ factored into left and right parts,

[TABLE]

As shown in Ref. [2] (Theorem 1) for boundary stochastic matrices, having nonvanishing both rates for stochastically setting the state of the site near the boundary, the propagator $U$ is irreducible and aperiodic. Hence, according to Perron-Frobenius theorem, the non-equilibrium steady state (NESS) eigenvector, corresponding to $\Lambda_{0}=1$, is unique and all other eigenvalues $\Lambda_{j\geq 1}$ are bounded by $|\Lambda_{j}|<1$ and thus the corresponding components of the state vector decay during the time evolution. These eigenvectors are also called decay modes, as they encode time evolution of any state as

[TABLE]

where $c_{j}$ are appropriate constants depending on the initial state.

In this paper we formulate a compact matrix product ansatz (MPA) which encodes eigenvectors of $U$ for the most important part of the spectrum, including NESS and the leading decay mode determining the spectral gap of $U$. In particular, we find that the spectrum organizes into orbitals, i.e., the sets of eigenvalues fulfilling the same Bethe equations (see Fig. 3 in Sec. 3.1). The NESS belongs to an especially simple (zeroth) orbital that contains also three additional eigenvectors whose eigenvalues do not depend on the system size. The first orbital contains the leading decay mode. The eigenvalues of the zeroth and the first orbital are nondegenerate. The other orbitals seem to be exponentially degenerate in system size.

We consider two types of boundary driving for which exact solutions can be found, the first we call conditional driving, and the second Bernoulli driving. Bernoulli driving has been introduced and studied for the steady state in Ref. [2].

1.1 Conditional boundary driving

For the conditional driving the boundary stochastic matrices read

[TABLE]

where $\alpha,\beta,\gamma,\delta\in[0,1]$ are some driving rates parametrizing the left and the right bath. We call this conditional driving since in $P^{\rm{L}}_{12}$ ($P^{\rm{R}}_{n-1,n}$) the probability of changing the site 1 depends only on the state of the neighboring site 2 (changing the site $n$ depends only on the state of the site $n-1$). For instance, if the site 2 is in state 0 then the site 1 will be stochastically set to state 0 with the rate $\alpha$ or to state 1 with the rate $1-\alpha$. If on the other hand, the site 2 is in state 1, the site 1 will be set to state 0 or 1 with the rates $\beta$ or $1-\beta$. The analogous also holds for $P^{\rm{R}}_{n-1,n}$.

1.2 Bernoulli boundary driving

For the Bernoulli driving, explained in more detail in Ref. [2], we have

[TABLE]

where $\alpha,\beta,\gamma,\delta\in[0,1]$ are again some driving rates specifying the left and the right bath. However, for this case it turns out that all the results are more compactly and conveniently expressed in terms of another set of, so called, difference parameters

[TABLE]

We note that both sets of boundary driving stochastic matrices (21,30) satisfy the commutativity condition (9) and represent so far the only known exactly solvable boundaries for the Rule 54. The converse seems not to be true. Not any pair of stochastic boundary matrices which satisfy (9) allow for exact solutions. Note as well that both type of boundary matrices (21,30) satisfy the conditions of ‘holographic ergodicity’ theorem of [2], implying exponential decay of any initial state to a unique NESS, for an open set of parameters $0<\alpha,\beta,\gamma,\delta<1$.

The rest of the paper is organized as follows. In Sec. 2 we introduce a cubic bulk algebra, which seems to provide the fundamental integrability relation of the model, and use it to solve the NESS-orbital in terms of MPA for both drivings (the NESS for the Bernoulli driving case was solved in terms of an alternative, patch-state ansatz in Ref. [2]). In Sec. 3 we introduce a generalization of the aforementioned algebra and use it to construct eigenvectors in the first orbital in terms of an inhomogeneous (spatially modulated) MPA. This orbital also contains the leading decay mode. The consistency conditions lead to a simple set of Bethe-like equations which yield the spectrum of the first orbital. We also discuss the thermodynamic limit and show that the spectral gap closes as $~{}1/n$. We close the section by discussing how the cubic bulk algebra may be re-written as a quadratic algebra, somewhat similar to Zamolodchikov-Faddeev (ZF) algebra. In Sec. 4 we provide a conjecture for the Bethe-like equations for the entire spectrum, as well as a conjecture for the degeneracy of the higher orbital eigenstates. Finally we end with the conclusions. The paper also contains an A stating explicitly all the components of the boundary vectors of the MPA generating the decay modes of the first orbital,.

Throughout the paper, whenever we provide explicit results, we will first state the results for the conditional driving (21) and then for the Bernoulli driving (30).

2 The cubic algebra and the non-equilibrium steady state

Let us first fix some notation. Quantities that are vectors in the physical space are written in bold-face. The numeral subscript of a physical space vector (or operator) denotes the site position in the tensor product $(\mathbb{R}^{2})^{\otimes n}$. When writing in component notation the components in physical space are labeled with a binary ‘spin’ index, such as $s\in\{0,1\}$. Matrices are written with capital roman letters. These are typically acting over 4-dimensional auxiliary space, except for the propagator $U$ and its local pieces $P_{k-1,k,k+1},P^{\rm{L}}_{12},P^{\rm{R}}_{n-1,n}$ which are operators in the physical space and act trivially in the auxiliary space. Matrices in an extended 8-dimensional auxiliary space (employed in the next section) will be denoted with hats. Row (column) vectors in the auxiliary space will be written as Dirac bras (kets).

We begin by defining a vector

[TABLE]

with components $W_{s}$ being $4\times 4$ matrices

[TABLE]

depending on a pair of formal parameters $\xi$ and $\omega$, which we call spectral parameters. We also define a matrix $\mathbf{W}^{\prime}$, which is $\mathbf{W}$ with $\xi$ and $\omega$ interchanged, i.e., (writing the dependence on the spectral parameters explicitly)

[TABLE]

The key property explored in this paper is a simple three-site cubic algebraic relation333Note a similar two-site cancellation mechanism in discrete time ASEP models [4, 5, 6]. which shall provide a cancellation mechanism to be used later for constructing the eigenstates of the Markov matrix

[TABLE]

where $S$ is a constant ‘delimiter’ matrix,

[TABLE]

satisfying $S^{2}=\mathbbm{1}_{4}$. By interchanging $\xi$ and $\omega$ in (34), i.e. interchanging $\mathbf{W}$ and $\mathbf{W}^{\prime}$, and multiplying with $P_{123}$ (noting that $P^{2}=\mathbbm{1}_{8}$) we obtain a dual bulk relation

[TABLE]

Note that Eq. (34) in fact represents $8$ matrix product identities, $W_{s}SW_{\chi(ss^{\prime}s^{\prime\prime})}W^{\prime}_{s^{\prime\prime}}=W_{s}W^{\prime}_{s^{\prime}}W_{s^{\prime\prime}}S$, and analogously for Eq. $(\ref{bulk1b})$, by writing out physical space components $s,s^{\prime},s^{\prime\prime}\in\{0,1\}$ for a vector on a triple of consecutive physical sites (denoted in (34,36) as $123$). Eqs. (34,36) thus represent a set of algebraic relations among $W_{0},W_{1},W^{\prime}_{0},W^{\prime}_{1},S$ which can be straightforwardly verified for the representation (32).

We shall begin by proposing a simple ansatz for the eigenvectors of $U$ in terms of the following staggered matrix product states

[TABLE]

In order for fixed point condition (11) to hold for $\mathbf{p},\mathbf{p^{\prime}}$ we require the following boundary conditions to be satisfied

[TABLE]

Specifically, writing out $U_{\rm e}\mathbf{p}$ in terms of (5) and the ansatz (37), one first uses Eq. (41) in order to introduce the delimiter $S$ in a string $\cdots\mathbf{W}^{\prime}_{n-3}\mathbf{W}_{n-2}\mathbf{W}^{\prime}_{n-1}S$ and then implements $\cdots P_{n-5,n-4,n-3}P_{n-3,n-2,n-1}$ via the dual bulk relation (36) in order to move the delimiter $S$ across to the left end where it is then absorbed, via $S^{2}=\mathbbm{1}$, by applying another boundary equation (40), arriving to $U_{\rm e}\mathbf{p}=\lambda_{\rm R}\mathbf{p}^{\prime}$. Analogously we proceed with $U_{\rm o}\mathbf{p}^{\prime}$, in terms of (6) and ansatz (38), now implementing the boundary Eqs. (39,42) to carry $S$ from left to right via the bulk relation (34), ending with $U_{\rm e}\mathbf{p}^{\prime}=\lambda_{\rm L}\mathbf{p}$. Thus, $\mathbf{p}$ is an eigenvector of $U$ with an eigenvalue444 Please note the use of small letters $\lambda,\lambda_{\rm L,R}$ for designating the spectral variables for the NESS-orbital in distinction to capitalised variables $\Lambda,\Lambda_{\rm L,R}$ referring to the general case (11) which, in the case of the first orbital, can be expressed as functions of the NESS-orbital data (see Sec. 3). $\lambda=\lambda_{\rm L}\lambda_{\rm R}$.

Solving the full set of boundary equations (39-42) should fix all the unknown parameters in the MPA (37,38) as the bulk relations are automatically satisfied. The solution is unique up to an irrelevant transformation of boundary vectors. We first focus on the conditional driving case (21) and then state the results for Bernoulli driving (30).

Solving separately the pair of boundary equations for the left side (39,40) we obtain the following unique solutions for the spectral parameters,

[TABLE]

and for the left boundary vectors (in physical space components)

[TABLE]

The right boundary equations (41,42), on the other hand, give the following unique solutions for the spectral parameters

[TABLE]

and the the right boundary vectors

[TABLE]

Now in order to get a consistent solution on both the left and right boundary we demand the two pairs of the spectral parameters $\xi$ and $\omega$ in Eqs. (43,44,46,47) to be equal. This gives us a closed pair of equations for $\lambda_{\rm L}$ and $\lambda_{\rm R}$,

[TABLE]

Rewriting these equations in terms of the eigenvalue $\lambda=\lambda_{\rm{L}}\lambda_{\rm{R}}$ leads to

[TABLE]

Clearly, $\lambda=1$ is always the solution, corresponding to NESS. The remainder is a cubic polynomial. Thus there are also three other solutions corresponding to three decay modes whose eigenvalues do not change with system size. This set of four eigenvalues shall be referred to as the NESS-orbital.

Following the same procedure for the case of Bernoulli driving (30) we find for the left boundary equations

[TABLE]

and on the right,

[TABLE]

yielding the characteristic polynomial for the eigenvalue

[TABLE]

The corresponding left boundary vectors are

[TABLE]

and the right boundary vectors are

[TABLE]

In summary, the NESS and three other decay modes whose eigenvalue does not depend on the system size can be obtained from the compatibility condition of the ‘scattering’ of MPA eigenvector (37,38) from the left and from the right stochastic boundary.

3 Generalization of the bulk algebra and the decay modes

The bulk algebra (34,36) admits several generalizations, one of which allows us to construct a set of decay modes for the two types of stochastic boundary drivings. For this purpose we introduce an additional parameter $z\in\mathbb{C}$, which will be referred to as momentum parameter, and define the following $4\times 4$ matrices

[TABLE]

We will also need define objects which are vectors in physical space and matrices in a 8-dimensional auxiliary space. We define block diagonal operators (with physical space component $s=0,1$),

[TABLE]

where $e_{ij}={|i\rangle}{\langle j|}$, $i,j\in\{1,2\}$, are $2\times 2$ unit matrices. Further define the following block triangular $8\times 8$ matrices in extended auxiliary space

[TABLE]

which depend on a pair of complex amplitude parameters $c_{+},c_{-}$. We now state the generalized inhomogeneous bulk relations

[TABLE]

where $\hat{S}=\mathbbm{1}_{2}\otimes S$, which can be straightforwardly checked to hold for any $\xi,\omega,z,c_{+},c_{-}\in\mathbb{C}$ and $k\in\mathbb{Z}$. Note that by setting $z=1,c_{+}=c_{-}=0$, we recover the original bulk algebra (34,36).

Defining a parity/swap transformation $\mathcal{R}:(\xi,\omega,z,c_{+},c_{-})\rightarrow(\omega,\xi,1/z,c_{-},c_{+})$, $\mathcal{R}^{2}=\mathrm{id}$, we find

[TABLE]

Applying $\mathcal{R}$ to the bulk relations (115) leads to another set of nonequivalent bulk relations (with $\hat{K}^{\prime(k)}:=\mathcal{R}\hat{K}^{(k)}$ and $\hat{L}^{\prime(k)}:=\mathcal{R}\hat{L}^{(k)}$). There are numerous other similar extensions of the bulk algebra (34,36), but they do not seem to be useful for constructing eigenvectors for the boundary drivings studied, so we omit writing them here.

Lemma: Let us assume that 8-dimensional boundary vectors ${\langle\hat{l}_{s}|}$, ${\langle\hat{l}^{\prime}_{s,s^{\prime}}|}$, ${|\hat{r}_{s,s^{\prime}}\rangle}$, ${|\hat{r}^{\prime}_{s}\rangle}$ exist, together with parameters $\xi,\omega,z,c_{+},c_{-},\Lambda_{\rm L},\Lambda_{\rm R}$, such that the following boundary equations are satisfied

[TABLE]

Then, the following inhomogeneous (site-dependent) MPA

[TABLE]

generates an eigenvector of $U=U_{\rm e}U_{\rm o}$ (5,6) with the eigenvalue $\Lambda=\Lambda_{\rm L}\Lambda_{\rm R}$.

Proof: Explaining how cancellation mechanism works is fully analogous to the simpler case of spatially homogeneous NESS-orbital (37-42). However, due to commutativity (9) we can explain it here for the reverse order: When $U_{\rm e}$ acts on $\mathbf{p}$, $P_{123}$ first acts on the left boundary vector and via (118) creates ${\hat{K}^{(1)}}{\hat{\mathbf{W}}_{3}}\hat{S}$. The subsequent $P$’s in $U_{\rm e}$ transfer ${\hat{K}^{(1)}}{\hat{\mathbf{W}}_{3}}\hat{S}$ to the right via the bulk algebra (115). Before the final $P_{n-3,n-2,n-1}$ acts, $P_{\rm{R}}$ acts (as it commutes with $P_{n-3,n-2,n-1}$) and creates the ${\hat{F}^{(n-3)}}{\hat{\mathbf{W}}^{\prime}_{1}}\hat{S}$ necessary for the final $P_{n-3,n-2,n-1}$ to transfer ${\hat{K}^{(n-5)}}{\hat{\mathbf{W}}_{n-3}}\hat{S}$ to the end as ${\hat{K}^{(n-3)}}{\hat{\mathbf{W}}_{n-1}}\hat{S}$, thus finally creating $\mathbf{p^{\prime}}$ (122). The odd-part of the propagator $U_{\rm o}$ acts analogously in reverse.

The MPA (121,122) will give us the leading decay mode and a set of other eigenvectors which we collectively call the first orbital. Due to the block upper triangular structure of $\hat{F}^{(k)},\hat{F}^{{}^{\prime}(k)},\hat{G}^{(k)},\hat{G}^{{}^{\prime}(k)}$, Eqs. (114), the matrix product state (121) (and analogously (122)) can be written as a superposition of terms containing a string of $W_{s}(\xi z,\omega/z)W^{\prime}_{s^{\prime}}(\xi z,\omega/z)$ and then $c_{\pm}z^{\pm k}G_{\pm}$ (or $c_{\pm}z^{\pm k}F_{\pm}$) at just before $W$ (or $W^{\prime}$) corresponding to site $k$ and then a string of $W_{s}(\xi/z,\omega z)W^{\prime}_{s^{\prime}}(\xi/z,\omega z)$. There is also a boundary term from $L_{\pm}$ in the superposition. Therefore, the first orbital can be understood as single quasiparticle excitations over the NESS, which are composed as superpositions of left- and right-propagating waves $z^{\pm k}$ with non-trivial scattering at the boundaries. This is somewhat similar in form to the matrix coordinate ansatz used in solving the decay modes of the ASEP model [7].

To solve the boundary equations we follow as similar procedure as outlined in Sec. 2 for the NESS-orbital. Let us decompose the extended auxiliary space as a direct sum of two 4-dimensional spaces ${\cal H}_{1}\oplus{\cal H}_{2}$, where an element of ${\cal H}_{i}$ is written as ${|i\rangle}\otimes{|\psi\rangle}$. Due to the upper triangular structure (114) the left boundary equations projected to the subspace ${\cal H}_{1}$ reduce to those for the NESS-orbital with a scaling factor $z$ coming from $e_{11}$ component of $\hat{L}^{(k)}$, but with rescaled spectral parameters ($\xi\to\xi z$, $\omega\to\omega/z$, because of (113)). Since the NESS solution we found in Sec. 2 is unique, the left boundary vector in this subspace must be the NESS-orbital boundary vector with scaled $\lambda_{\rm L}=\Lambda_{\rm L}/z$. Comparing (117,118) with (39,40) and (43,44) immediately fixes the values of the spectral parameters $\xi$ and $\omega$ (writing first for the conditional driving (21)):

[TABLE]

The remaining components of the left boundary equations in the subspace ${\cal H}_{2}$ come from either the diagonal components $e_{22}$ or the off-diagonal components $e_{12}$ of the auxiliary space operators. Requiring that the equations are solved for arbitrary $z,\alpha,\beta$, this fixes the ratio of the amplitudes

[TABLE]

and that the general form of the left boundary vectors must be

[TABLE]

where the explicit form of the ‘offdiagonal’ vectors ${\langle\tilde{l}_{s}|}$ and ${\langle\tilde{l}^{\prime}_{s,s^{\prime}}|}$ are given in A. Following a fully analogous procedure for the right boundary equations (119), (120), we arrive to the following results for the spectral parameters and the amplitude ratio

[TABLE]

where $m=\frac{n}{2}-2$, and for the right boundary vectors

[TABLE]

whose components in ${\cal H}_{2}$ are expressed in terms of right boundary vectors for the NESS-orbital and the complementary (offdiagonal) components ${|\tilde{r}_{s,s^{\prime}}\rangle}$ and ${|\tilde{r}^{\prime}_{s}\rangle}$ are given in A. Pairwise identifying equations (123,127), (124,128), and (125,129) we obtain a closed set of Bethe-like equations for $\Lambda_{\rm L},\Lambda_{\rm R},z$:

[TABLE]

Writing $\Lambda_{\rm L}=\Lambda/\Lambda_{\rm R}$ and eliminating the variable $\Lambda_{\rm R}$ these equations can be rewritten as a pair of algebraic equations for $\Lambda,z$, the first of which can be understood as a nonequilibrium quasiparticle dispersion relation

[TABLE]

This equation has $4(2m+1)$ distinct roots which comprise the first orbital.

For the case of Bernoulli driving (30) we find the following solutions for the left boundary equations (117,118)

[TABLE]

and for the right ones (119,120)

[TABLE]

Note that in this case these equations only depend on the difference of the driving parameters $\mu=\alpha-\beta$ (on the left) and $\sigma=\gamma-\delta$ (on the right) and thus so will the eigenvalues. These lead to the following set of of Bethe-like equations,

[TABLE]

or equivalently

[TABLE]

The boundary vectors in the Bernoulli driving case are of the same form as in the conditional driving case, namely (126,130). The explicit expressions for their components are given in A.

3.1 Thermodynamics

In this subsection we will study the thermodynamic properties of the eigenvalues of the first orbital. Expanding equation (135) (or (146)) in $1/m$ (recalling that $m=\frac{n}{2}-2$) gives us that in the leading order of $1/n$ ($(1/n)^{0}=1)$, $z$ must be unimodular

[TABLE]

where $\kappa$ is the quasi-momentum which may be restricted to interval $[0,\pi)$ as (134) and (135) (or (145) and (146)) are symmetric under the transformation $z\to-z$. Thus, in the leading order of $1/n$ expansion the spectrum in the first orbital converges to an algebraic curve $\Lambda(\kappa)$ given by (147) and (134) (or (145)) (see Fig. 3)

The case $\Lambda(\kappa=0)=1$ corresponds to the eigenvalues which in the thermodynamic limit converge to $1$ and their scaling with $1/n$ will determine the asymptotic relaxation rate of the system. Writing the expansion of $\Lambda(0)$ as

[TABLE]

and inserting into (134) (and (145)) we find that the gap closes as $1/n$ for both drivings. In the case of conditional driving (21) we find

[TABLE]

and in the case of Bernoulli driving (30) we find

[TABLE]

where $\mu=\alpha-\beta$ and $\sigma=\gamma-\delta$. This $1/n$ scaling of the gap is consistent with the results of Ref. [2] and ballistic transport.

3.2 Quadratic form of the bulk algebra

Even though the most elementary, partonic blocks of our exact solutions satisfy cubic algebraic relations, either (34,36) or (115), it is possible to rewrite them in terms of a quadratic algebra at the expense of defining auxiliary space operators which depend on two adjacent physical sites rather than one.

Let us define the following inhomogeneous 4-component vectors of $8\times 8$ matrices

[TABLE]

The cubic bulk algebra (115) is then equivalent to the following quadratic algebra (formulated as 16-component vectors over four consecutive physical sites $1234$)

[TABLE]

This quadratic formulation has perhaps some appeal as it is reminiscent of ZF algebra.555Essential differences to any meaningful formulation of ZF algebra still remain, most notably, our scattering operator $P$ has no ‘momentum’ dependence. Perhaps this can be mended by some stochastic deformation of the model. The eigenvector MPA (121,122) now transforms to:

[TABLE]

where, as always $m=n/2-2$ and,

[TABLE]

Note that one can formulate and solve the boundary equations solely in terms of such two-site boundary vectors arriving to equivalent Bethe-like equations as in the previous section.

4 Conjectures about the complete spectrum: modified Bethe-like equations and the degeneracy

Despite numerous attempts we were not able to construct any other eigenvectors of $U$ beyond NESS-orbital and the first orbital. Still, numerical inspections of the spectrum (see Fig. 3) suggest that the structure of higher orbitals should be very similar, but the eigenvalues become degenerate with the degeneracy which quickly increases with the level of the orbital. These comprise in total $2^{n-2}$ nonvanishing eigenvalues while the eigenvalue 0 is $3\times 2^{n-2}$ fold degenerate.

We were able to guess the Bethe-like equations which reproduce the entire spectrum. Introducing an integer $p$ which counts the orbital level and runs from $1$ to $m=n/2-2$, we postulate the modified Bethe-like equations, either for the conditional driving (21):

[TABLE]

or for the Bernouli driving (30):

[TABLE]

Each of these orbitals has exactly $4(2m+1)$ distinct roots. We conjecture that the above equations, together with the NESS characteristic polynomial (51,56), describe the entire spectrum of $U$. Indeed, agreement with the spectrum obtained from numerical diagonalization of $U$, for up to $n=16$, is perfect, within trustable precision of numerical routines (as demonstrated in Fig. 3).

Furthermore, we provide a conjecture for the average degeneracy of the eigenvalues of the $p$-th orbital which seems to scale as

[TABLE]

Note that the function $g(m,p)$ can also be non-integer for some values of $m,p$. This means that there are different degeneracies within the same orbital $p$ for system size $n=2(m+2)$ and $g(m,p)$ gives the average value of the degeneracy. This conjecture has also been confirmed numerically for up to $n=16$. We can make a simple consistency check by counting the total number of roots in all the orbitals together with their multiplicity $\sum_{p=1}^{m}4(2m+1)g(m,p)=4(4^{m}-1)=2^{n-2}-4$ which together with four eigenvalues of the NESS-orbital yield the total number of $2^{n-2}$ nonvanishing eigenvalues.

5 Conclusions and open problems

We found that the spectrum of the Markov matrix of a deterministic boundary driven cellular automaton (the Rule 54) organizes into orbitals. Throughout the paper we gave results for two types of non-equivalent stochastic boundary drivings, (21) and (30). The eigenvalues in each orbital fulfil a set of three coupled Bethe-like equations. We found explicit matrix product forms of the eigenvectors in two main orbitals – the NESS-orbital (containing the NESS and three other eigenvectors) and the first orbital (which contains the leading decay mode and eigenvalues of the largest modulus). To find the NESS-orbital we used a 4-dimensional representation (with two spectral parameters) of a three-site bulk algebra cancellation mechanism (34). This three-site bulk algebra is similar in form to a two-site bulk cancellation mechanism in discrete time ASEP models [4, 5, 6]. To find the first orbital we generalized the aforementioned bulk algebra to an 8-dimensional auxiliary space (115). The structure of such positionally dependent bulk algebra allows for construction of the eigenvectors in a compact form of an inhomogeneous MPA, which may be understood as a superposition of local single-particle excitations over the NESS-orbital with different momenta, reminiscent of the matrix coordinate ansatz for ASEP models [7]. We also investigated the thermodynamic properties and proved that the spectral gap yielding the ultimate relaxation time scales as $~{}1/n$. In the thermodynamic limit, the first orbital defines an algebraic curve in the complex plane which borders the spectrum of the model. We have also shown how the cubic bulk algebra may be rewritten as a quadratic algebra with operators acting on two physical sites, instead of one. We note that our inhomogeneous MPA can be rewritten as well in the form of an inhomogeneous patch state ansatz as proposed and used to find NESS of the model in Ref. [2], however the elements of the patch tensors have to be replaced by positionally dependent $2\times 2$ matrices. This has been actually the route through which we arrived at the results reported in the present paper.

Although we were unable to find all the eigenvectors, we provided a conjecture for the Bethe equations for all the orbitals, which has been corroborated by extensive numerical investigation (Sec. 4). The higher orbitals are highly degenerate and we provide a conjecture for the average degeneracy of each orbital in the same section.

Many open questions remain. Most pressingly, one would wish to construct the eigenvectors in all the orbitals exactly and prove the conjecture for the general Bethe-like equations as given in Sec. 4. A direct generalization of the bulk algebra (115) to two-particle excitations does not seem to work, and numerous other generalizations are possible making it difficult to ascertain how to continue. The fact that the higher orbital Bethe-like equations are very similar to the first orbital’s might suggest that the nonequilibrium quasiparticle excitations are non-interacting and that $z$ represents simple the center-of-mass momentum of all excitations. However, such an idea seems to be an oversimplification since the independent particle model cannot explain exponentially large (in orbital level) degeneracy of the eigenvalues. We thus interpret this as an interacting model with an exponential bunching of quasiparticles. It is also not clear at present if and how the exact solution of the boundary driven Rule 54 model reported here connects to Yang-Baxter equation and common language of integrability.

Another interesting question which remains is whether the cases (21) and (30) complete the set of integrable stochastic boundaries of the model or not. The problem of classification of integrable stochastic boundaries of a bulk deterministic (e.g. Hamiltonian) integrable theory is generally open.

Even though we were unable to construct the complete set of eigenvectors of our model we feel that the results obtained here can be extended for use in other models, for instance, to complement the study of asymptotic decay of densities in reaction models [8], help in the study of the decay modes of discrete time models, notably discrete time ASEP models [4, 5, 6], and in particular, to other driven integrable cellular automata with deterministic bulk dynamics (e.g., [3, 9, 10, 11, 12]) and perhaps even driven quantum models in discrete time [13].

Acknowledgements

We thank V. Popkov, M. Vanicat, E. Ragoucy for discussions and E. Ilievski and L. Zadnik for useful comments on the manuscript. The work has been supported by the grants P1-0044 and N1-0025 of Slovenian Research Agency, and ERC Grant OMNES.