Integrability of a deterministic cellular automaton driven by stochastic boundaries

In 1993, Bobenko et al [1] presented and discussed a simple model of a two-state (${{\mathbb{Z}}}_{2}$) reversible cellular automaton—the so-called 'rule 54'—(RCA54) which possesses the main features of integrability and can perhaps be considered as the simplest strongly interacting classical dynamical system with nontrivially scattering soliton solutions. The rule can be encoded via 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}$, determining the state of the south edge in terms of the states of the west, east, and north edges

$$\begin{eqnarray}&&{s}_{{\rm{S}}}=\chi ({s}_{{\rm{W}}},{s}_{{\rm{N}}},{s}_{{\rm{E}}})={s}_{{\rm{N}}}+{s}_{{\rm{W}}}+{s}_{{\rm{E}}}+{s}_{{\rm{W}}}{s}_{{\rm{E}}}\;(\mathrm{mod}\;\;2),\quad {s}_{{\rm{S}}},{s}_{{\rm{N}}},{s}_{{\rm{W}}},{s}_{{\rm{E}}}\in {{\mathbb{Z}}}_{2},\end{eqnarray} \tag{ 1 }$$

where the time runs in north–south direction.

The RCA54 is clearly reversible, as at the same time ${s}_{{\rm{N}}}=\chi ({s}_{{\rm{W}}},{s}_{{\rm{E}}},{s}_{{\rm{S}}})$. See the bulk (central) part of figure 1 to observe a typical evolution pattern of the RCA54 rule when applied to some initial configuration in the upper rows. The code of such cellular automaton is given by an ordered list of binary digits $(\chi (000),\chi (001)$, $\chi (010),\chi (011)$, $\chi (100),\chi (101)$, $\chi (110)$, $\chi (111))$ for all possible 2³ values of binary triples $({s}_{{\rm{W}}},{s}_{{\rm{N}}},{s}_{{\rm{E}}})$, namely ${\sum }_{s,s^{\prime} ,s^{\prime\prime} }\chi (s,s^{\prime} ,s^{\prime\prime} ){2}^{4s+2s^{\prime} +s^{\prime\prime} }=54$.

Zoom In Zoom Out Reset image size

Deterministic long-time dynamics generated by (1) via permutations over the set ${ \mathcal C }={({{\mathbb{Z}}}_{2})}^{\times n}$ of all possible ${2}^{n}$ configurations of the zig-zag lattice (chain) of n subsequent cells (see e.g. figure 2 for schematic illustration) is always, either rather trivial, or depending crucially on the boundary conditions, i.e. the states of the cells at spatial coordinates x = 1 and x = n which cannot be determined dynamically unless additional rules are introduced. In this paper we propose stochastic update rules for the boundary cells implemented through a pair of local Markov chains which can be interpreted as chemical reservoirs parametrized with absorbtion and emission rates of solitons at each boundary. Note that such a hybrid bulk-deterministic, boundary-stochastic statistical mechanics paradigm is fundamentally different from well-studied boundary driven-diffusive classical simple exclusion processes [2], which have stochastic bulk dynamics but can be interpreted as many-body Markov chains over an identical state space ${{\mathbb{R}}}^{{ \mathcal C }}\simeq {({{\mathbb{R}}}^{2})}^{\otimes n}$. Rather, the novel paradigm proposed here should be understood as the simplest classical-dynamical version of integrable boundary driven nonequilibrium quantum spin chains which have recently been intensely studied (see e.g. [9] for a review).

Zoom In Zoom Out Reset image size

After introducing the Markov chain model for boundary driven RCA54 dynamics in section 2, we shall present the proof of irreducibility and aperiodicity of the resulting Markov matrix, which implies uniqueness of the nonequilibrium steady state (NESS) and asymptotic approach to NESS from any initial state (ergodicity and dynamical mixing—sometimes referred to as strong ergodicity [7]). The main result of our paper, presented in section 3, is an exact and explicit solution for NESS in terms of a particular, commutative but correlated (nonseparable) matrix product ansatz, which we term a patch state ansatz (PSA). In subsequent section 4 we shall demonstrate explicit computation of observables, such as steady state values of density, soliton-currents, and arbitrary k-point spatial density–density correlation functions. Moreover, we demonstrate in section 5, that similarly to the case of boundary driven quantum XXZ spin-1/2 chains [8, 9], one can exploit the analytical form of NESS to generate nontrivial (quasi)local conservation laws. In the last section 6 we discuss some interesting follow-up questions, such as computation of decay modes and full spectrum of the Markov matrix, and conclude.

Throughout our work we shall—for simplicity and symmetry reasons—assume that the number of cells is even

$$\begin{eqnarray}&&n=2m.\end{eqnarray} \tag{ 2 }$$

The extension of the results to the case of odd sizes n should be straightforward. We identify a state space with a vector space ${ \mathcal S }={{\mathbb{R}}}^{{ \mathcal C }}={({{\mathbb{R}}}^{2})}^{\otimes n}$, a linear space embedding a convex subspace of all probability state vectors $\underline{p}=({p}_{0},{p}_{1},\ldots ,{p}_{{2}^{n}-1})\in { \mathcal S }$, satisfying ${p}_{s}\geqslant 0,{\sum }_{s=0}^{{2}^{n}-1}{p}_{s}=1$. A vector of k binary digits $\underline{s}=({s}_{1},{s}_{2},\ldots ,{s}_{k})\in {{\mathbb{Z}}}_{2}^{\times k}$ shall often be identified with an integer $s={\sum }_{j=1}^{k}{s}_{k}{2}^{k-j}$, or, components of the probability state vector shall be written as ${p}_{s}\equiv {p}_{({s}_{1},{s}_{2},\ldots ,{s}_{n})}\equiv {p}_{{s}_{1},{s}_{2},\ldots {s}_{n}}$. Deterministic, local RCA54 rule in the bulk (1) can be encoded into a ${2}^{3}\times {2}^{3}$ permutation matrix

$$\begin{eqnarray}&&{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} },\end{eqnarray} \tag{ 3 }$$

or

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

which is self invertible, ${P}^{2}={{\mathbb{1}}}_{{2}^{3}}$. Here, ${{\mathbb{1}}}_{d}$ denotes a d-dimensional identity matrix and ${\delta }_{s,t}$ a Kronecker symbol.

On the boundaries, we define two-site local stochastic Markov chains, which depend on the state of a pair of near boundary cells. Firstly, we define a simple, single-cell (ultralocal) Markov chain, depending on the conditional probability α to change cell $0\to 1$ and conditional probability β to change cell $1\to 0$

$$\begin{eqnarray}{E}^{\alpha ,\beta }=\left(\begin{array}{cc}1-\alpha & \beta \\ \alpha & 1-\beta \end{array}\right),\quad \alpha ,\beta \in [0,1].\end{eqnarray} \tag{ 4 }$$

Secondly, we define two-cell local Markov chains for a pair of cells near the boundary by imagining another, stochastic cell just beyond the boundary following a Bernoulli process $B\left(\frac{1}{2},\frac{1}{2}\right)$ and then applying the local RCA54 rule (1) to the triple of cells. Such processes are generated by the following 4 × 4 Markov matrices matrices, for each boundary

$$\begin{eqnarray}{\tilde{P}}_{(s^{\prime} ,s^{\prime\prime} ),(t^{\prime} ,t^{\prime\prime} )}^{{\rm{L}}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{s=0}^{1}{P}_{(s,s^{\prime} ,s^{\prime\prime} ),(s,t^{\prime} ,t^{\prime\prime} )},\\ {\tilde{P}}_{(s,s^{\prime} ),(t,t^{\prime} )}^{{\rm{R}}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{s^{\prime\prime} =0}^{1}{P}_{(s,s^{\prime} ,s^{\prime\prime} ),(t,t^{\prime} ,s^{\prime\prime} )}.\end{eqnarray} \tag{ 5 }$$

These relations can be compactly written in terms of a partial trace ${{\rm{tr}}}_{k}$ over kth qubit of ${({{\mathbb{C}}}^{2})}^{\otimes 3}$,

$$\begin{eqnarray*}&&{\tilde{P}}^{{\rm{L}}}=\displaystyle \frac{1}{2}{{\rm{tr}}}_{1}P,\quad {\tilde{P}}^{{\rm{R}}}=\displaystyle \frac{1}{2}{{\rm{tr}}}_{3}P.\end{eqnarray*}$$

Composing these two Markov processes, we obtain, for each boundary, the final forms of 4 × 4 matrices of two-cell boundary Markov chains

$$\begin{eqnarray}{P}^{{\rm{L}}} & = & {\tilde{P}}^{{\rm{L}}}({E}^{\alpha ,\beta }\otimes {{\mathbb{1}}}_{2})=\left(\begin{array}{cccc}\frac{1}{2} & & \frac{1}{2} & \\ & 1-\alpha & & \beta \\ \frac{1}{2} & & \frac{1}{2} & \\ & \alpha & & 1-\beta \end{array}\right),\\ {P}^{{\rm{R}}} & = & {\tilde{P}}^{{\rm{R}}}({{\mathbb{1}}}_{2}\otimes {E}^{\gamma ,\delta })=\left(\begin{array}{cccc}\frac{1}{2} & \frac{1}{2} & & \\ \frac{1}{2} & \frac{1}{2} & & \\ & & 1-\gamma & \delta \\ & & \gamma & 1-\delta \end{array}\right).\end{eqnarray} \tag{ 6 }$$

The full Markov chain propagator $U\in {\rm{End}}({ \mathcal S })$ is then written as a composition of two temporal layer propagators

$$\begin{eqnarray}&&U={U}_{{\rm{o}}}{U}_{{\rm{e}}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{U}_{{\rm{e}}}={P}_{123}{P}_{345}\cdots {P}_{n-3,n-2,n-1}{P}_{n-1,n}^{{\rm{R}}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{U}_{{\rm{o}}}={P}_{n-2,n-1,n}\cdots {P}_{456}{P}_{234}{P}_{12}^{{\rm{L}}},\end{eqnarray} \tag{ 9 }$$

where embedding into ${\rm{End}}({ \mathcal S })$ is understood as ${P}_{k,k+1,k+2}={{\mathbb{1}}}_{{2}^{k-1}}\otimes P\otimes {{\mathbb{1}}}_{{2}^{n-k-2}}$, ${P}_{1,2}^{{\rm{L}}}={P}^{{\rm{L}}}\otimes {{\mathbb{1}}}_{{2}^{n-2}},{P}_{n-1,n}^{{\rm{R}}}={{\mathbb{1}}}_{{2}^{n-2}}\otimes {P}^{{\rm{R}}}$. The model depends on four external driving parameters $\alpha ,\beta ,\gamma ,\delta \in [0,1]$, which can be understood (or related to) injection/absorption rates of solitons at the left/right boundary, respectively. Monte Carlo dynamics of driven RCA54 for some typical values of driving parameters is illustrated in figure 1, while schematic composition of the full many-body Markov generator is depicted in figure 2.

The many-body propagator U is clearly a stochastic matrix, i.e., its non-negative elements in each column sum to 1. In fact, for generic values of driving parameters $0\lt \alpha ,\beta ,\gamma ,\delta \lt 1$, it has exactly 4 nonvanishing matrix elements in each column, corresponding to four combinations of boundary cells x = 1 and x = n. Our goal is to find a steady state probability state vector $\underline{p}\in { \mathcal S }$, which is nothing but a fixed point of our many-body Markov chain—the NESS:

$$\begin{eqnarray}&&U\underline{p}=\underline{p},\end{eqnarray} \tag{ 10 }$$

or equivalently, to find a pair of probability state vectors $\underline{p},\underline{p}^{\prime} \in { \mathcal S }$ on subsequent temporal zig-zag layers, satisfying

$$\begin{eqnarray}&&{U}_{{\rm{e}}}\underline{p}=\underline{p^{\prime} },\quad {U}_{{\rm{o}}}\underline{p^{\prime} }=\underline{p}.\end{eqnarray} \tag{ 11 }$$

However, establishing an existence of a unique NESS and relaxation towards NESS from an arbitrary initial probability state vector amounts to [7] showing the following statement:

Theorem 1. The ${2}^{n}\times {2}^{n}$ matrix U, equation (7), is irreducible and aperiodic for generic values of driving parameters, more precisely, for an open set $0\lt \alpha ,\beta ,\gamma ,\delta \lt 1$.

Proof. We recall [5] that a finite, non-negative matrix U is irreducible if for any pair of configurations $\underline{s},\underline{s}^{\prime} \in { \mathcal C }$, one can find a natural number ${t}_{0}\in {\mathbb{N}}$ such that ${({U}^{{t}_{0}})}_{\underline{s}^{\prime} ,\underline{s}}\gt 0$. An irreducible matrix U is aperiodic if for some configuration $\underline{s}\in { \mathcal C }$, the greatest common divisor of recurrence times $t\in {\mathbb{N}}$, for which ${({U}^{t})}_{\underline{s},\underline{s}}\gt 0$, is 1.

Let us first show aperiodicity. As we have argued above, ${U}_{\underline{s}^{\prime} ,\underline{s}}$ connects each configuration $\underline{s}$ with exactly 4 other configurations $\underline{s}^{\prime}$ with all possible values of boundary cells, $({s}_{1}^{\prime },{s}_{n}^{\prime })\in \{(0,0),(0,1),(1,0),(1,1)\}$, unless some of the parameters $\alpha ,\beta ,\gamma ,\delta$ is equal exactly 0 or 1 which is the marginal case that is excluded from the discussion. Let us now take a sufficiently large positive integer t₀, to be determined below, and fix $\underline{s}({t}_{0})\equiv \underline{s}^{\prime} ,\underline{s}(0)\equiv \underline{s}$. We shall then construct a walk

$$\begin{eqnarray}&&\underline{s}(0)\to \underline{s}(1)\to \underline{s}(2)\to \cdots \;\to \;\underline{s}({t}_{0}),\end{eqnarray} \tag{ 12 }$$

i.e., a path through the Markov graph defined by positive elements of U, which connects $\underline{s}$ and $\underline{s}^{\prime}$ in t₀ steps and implies ${({U}^{{t}_{0}})}_{\underline{s}^{\prime} ,\underline{s}}\gt 0$ (see figure 3 for a 'self-contained' graphic illustration of the idea of proof). We are still free to choose the values of the boundary cells ${s}_{1,n}(t)$ along the walk $t\in \{1,2,\ldots ,{t}_{0}-1\}$ apart from the ends. For the first part of the walk $t=1,2\ldots {t}_{+}$, up to some ${t}_{+}\lt {t}_{0}$, we are fixing them with the rule

$$\begin{eqnarray}&&{s}_{1}(t)={s}_{2}(t-1),\quad {s}_{n}(t)={s}_{n-1}(t-1).\end{eqnarray} \tag{ 13 }$$

The evolution of the interior values of the cells s_x(t) for $1\lt x\lt n$ and $t\leqslant {t}_{+}$ is then completely specified by the deterministic RCA54, while (13) provide the causal absorbing boundary conditions. Indeed, each time the boundary cell, say x = 1, gets occupied, ${s}_{1}(t)=1$, the soliton as defined in [1] is absorbed (see figure 3). As the solitons only move ballistically (at speed 1) and scatter pairwise (with time-lag 1), while they cannot form bound states, it is clear that a finite time scale ${t}_{+}\in {\mathbb{N}}$ exists, surely smaller than n², after which all the solitons will be absorbed and we end up in a vacuum configuration $\underline{s}({t}_{+})=(0,0\ldots ,0)$.

For the rest of the walk $t\in \{{t}_{+}+1,\ldots {t}_{0}\}$ we need to show that an alternative boundary rules exist which create the configuration $\underline{s}^{\prime}$ out of the vacuum in another ${t}_{-}={t}_{0}-{t}_{+}$ steps. This is easily achieved by using time-reversibility of RCA54 and arguing that a vacuum configuration is again generated from $\underline{s}^{\prime}$ in some ${t}_{-}$ steps if the anti-causal absorbing boundary conditions are set (which are equivalent to (13) when the time runs backwards)

$$\begin{eqnarray}&&{s}_{1}(t)={s}_{2}(t+1),\quad {s}_{n}(t)={s}_{n-1}(t+1),\quad {\rm{for}}\quad t={t}_{0}-1,{t}_{0}-2,\ldots ,{t}_{0}-{t}_{-}.\end{eqnarray} \tag{ 14 }$$

The entire walk then connects $\underline{s}$ to $\underline{s}^{\prime}$ in ${t}_{0}={t}_{+}+{t}_{-}$ steps and implies ${({U}^{{t}_{0}})}_{\underline{s}^{\prime} ,\underline{s}}\gt 0$, for arbitrary pair $\underline{s},\underline{s}^{\prime} \in { \mathcal C }$ where the minimal possible integer t₀ may depend on the choice of $\underline{s},\underline{s}^{\prime}$. This proves irreducibility of (7).

Considering $\underline{s}^{\prime} =\underline{s}$, we have just shown that ${U}_{\underline{s},\underline{s}}^{{t}_{0}}\gt 0$ for some t₀ depending on $\underline{s}$. But since in between annihilating the configuration $\underline{s}$ in ${t}_{+}$ time steps and then creating it again in another ${t}_{-}$ steps³ , while ${t}_{0}={t}_{+}+{t}_{-}$, we can await in the vacuum state for an arbitrary additional number of steps ${t}_{{\rm{gap}}}\geqslant 0$, i.e. increase the walk by a segment of ${t}_{{\rm{gap}}}$ intermediate vacuum configurations, and still have ${({U}^{{t}_{0}+{t}_{{\rm{gap}}}})}_{\underline{s},\underline{s}}\gt 0$. The greatest common divisor of the set $\{{t}_{0}+{t}_{{\rm{gap}}};{t}_{{\rm{gap}}}\in {{\mathbb{Z}}}_{+}\}$ is clearly 1, so we have shown aperiodicity.□

In fact, a careful combinatorics of soliton scatterings and boundary absorbtions reveals that the minimal time t₀ which suffices for all pairs of configurations $\underline{s},\underline{s}^{\prime}$, i.e. after which ${{\rm{U}}}^{{t}_{0}}$ becomes a (strictly) positive matrix, reads

$$\begin{eqnarray}&&\mathrm{min}\{{t}_{0}\in {\mathbb{N}};{({U}^{{t}_{0}})}_{\underline{s},\underline{s}^{\prime} }\gt 0,\forall \underline{s},\underline{s}^{\prime} \}=\displaystyle \frac{3}{2}n-2.\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

In conclusion, the Perron–Frobenius theorem [5] guaranties that NESS probability state vector $\underline{p}$ satisfying the fixed point condition (10), or (11), is unique—eigenvalue 1 of U is simple—and all other eigenvalues of U lie strictly inside the unit circle. As a consequence, the Markov dynamics $\underline{p}(t)={U}^{t}\underline{p}(0)$ is ergodic and mixing and an arbitrary initial probability state vector $\underline{p}(0)$ converges exponentially in t to NESS.

We shall now explicitly construct the probability state vectors $\underline{p}$ and $\underline{p}^{\prime}$ of NESS, solving equation (11) in terms of a simple ansatz, which we term a PSA (illustrated in figure 4).

Zoom In Zoom Out Reset image size

Theorem 2. For an open set of driving parameters, $0\lt \alpha ,\beta ,\gamma ,\delta \lt 1$, the NESS solution $\underline{p},\underline{p}^{\prime} \in { \mathcal S }$ of the fixed point condition (11) can be written, for any even size n, in the form

$$\begin{eqnarray}{p}_{{s}_{1},{s}_{2},\ldots ,{s}_{n}} & = & {L}_{{s}_{1}{s}_{2}{s}_{3}}{X}_{{s}_{2}{s}_{3}{s}_{4}{s}_{5}}{X}_{{s}_{4}{s}_{5}{s}_{6}{s}_{7}}\cdots {X}_{{s}_{n-4}{s}_{n-3}{s}_{n-2}{s}_{n-1}}{R}_{{s}_{n-2}{s}_{n-1}{s}_{n}},\\ {p}_{{s}_{1},{s}_{2},\ldots ,{s}_{n}}^{\prime } & = & {L}_{{s}_{1}{s}_{2}{s}_{3}}^{\prime }{X}_{{s}_{2}{s}_{3}{s}_{4}{s}_{5}}^{\prime }{X}_{{s}_{4}{s}_{5}{s}_{6}{s}_{7}}^{\prime }\cdots {X}_{{s}_{n-4}{s}_{n-3}{s}_{n-2}{s}_{n-1}}^{\prime }{R}_{{s}_{n-2}{s}_{n-1}{s}_{n}}^{\prime },\end{eqnarray} \tag{ 16 }$$

for some rank-4 and rank-3 tensors of strictly positive components ${X}_{{ss}^{\prime} {uu}^{\prime} },{X}_{{ss}^{\prime} {uu}^{\prime} }^{\prime },{L}_{{suu}^{\prime} },{L}_{{suu}^{\prime} }^{\prime },{R}_{{ss}^{\prime} u},{R}_{{ss}^{\prime} u}^{\prime }$, with binary indices $s,s^{\prime} ,u,u^{\prime} \in \{0,1\}$.

Explicit n-independent algebraic expressions for the tensors $X,X^{\prime} ,L,L^{\prime} ,R,R^{\prime}$ in terms of the parameters of the model $\alpha ,\beta ,\gamma ,\delta$ shall be given later in the proof (20) and (23).

Proof. We shall first (i) present a minimal set of equations which are sufficient to determine the tensors $X,X^{\prime} ,L,L^{\prime} ,R,R^{\prime}$ under the assumption of the theorem. These nonlinear algebraic equations, being sufficiently simple, can be readily solved. Then, in the second part of the proof (ii) we shall show that the PSA solution identically satisfies every component of the fixed point conditions (11), for any even n.

(i) A normalization of the PSA (16) can be chosen such that

$$\begin{eqnarray}&&{X}_{0000}={X}_{0000}^{\prime }=1,\quad {L}_{000}={R}_{000}^{\prime }=1.\end{eqnarray} \tag{ 17 }$$

Clearly ${X}_{0000}={X}_{0000}^{\prime }$, otherwise the probabilities of the vacuum configurations ${p}_{\mathrm{0,0},\ldots 0}$ and ${p}_{0,0,\ldots ,0}^{\prime }$ would scale differently with n which is not possible since ${U}_{{\rm{o}}}$ and ${U}_{{\rm{e}}}$ directly connect $(0,0,\ldots ,0,0)$ only with configurations $({s}_{1},0,\ldots ,0,{s}_{n})$.

Let us now assume the ansatz (16) and write all components of equations (11), ${({U}_{{\rm{e}}}\underline{p}-\underline{p}^{\prime} )}_{\underline{s}}={({U}_{{\rm{o}}}\underline{p}^{\prime} -\underline{p})}_{\underline{s}}=0$, pertaining to four-cluster configurations in the bulk of the form⁴  $\underline{s}=({0}^{\{2k+1\}},s,s^{\prime} ,u,u^{\prime} ,{0}^{\{n-5-2k\}})$, for $k=0,1,\ldots ,m-3$, and three-cluster configurations near each boundary, $\underline{s}=(v^{\prime} ,s,s^{\prime} ,{0}^{\{n-3\}})$ and $\underline{s}=({0}^{\{n-3\}},s,s^{\prime} ,u)$, resulting in the following finite set of equations⁵ :

$$\begin{eqnarray}&&{X}_{00{ss}^{\prime} }^{\prime }{X}_{{ss}^{\prime} {uu}^{\prime} }^{\prime }{X}_{{uu}^{\prime} 00}^{\prime }{X}_{0000}^{\prime }\\ &&\quad ={X}_{00\chi (0{ss}^{\prime} )s^{\prime} }{X}_{\chi (0{ss}^{\prime} )s^{\prime} \chi (s^{\prime} {uu}^{\prime} )u^{\prime} }{X}_{\chi (s^{\prime} {uu}^{\prime} )u^{\prime} \chi (u^{\prime} 00)0}{X}_{\chi (u^{\prime} 00)000},\\ &&{X}_{0000}{X}_{00{ss}^{\prime} }{X}_{{ss}^{\prime} {uu}^{\prime} }{X}_{{uu}^{\prime} 00}\\ &&\quad ={X}_{000\chi (00s)}^{\prime }{X}_{0\chi (00s)s\chi ({ss}^{\prime} u)}^{\prime }{X}_{s\chi ({ss}^{\prime} u)u\chi ({uu}^{\prime} 0)}^{\prime }{X}_{u\chi ({uu}^{\prime} 0)00}^{\prime },\\ &&{L}_{v^{\prime} {ss}^{\prime} }^{\prime }{X}_{{ss}^{\prime} 00}^{\prime }{X}_{0000}^{\prime }{R}_{000}^{\prime }={L}_{v^{\prime} \chi (s^{\prime} {ss}^{\prime} )s^{\prime} }{X}_{\chi (v^{\prime} {ss}^{\prime} )s^{\prime} \chi (s^{\prime} 00)0}{X}_{\chi (s^{\prime} 00)000}\displaystyle \frac{{R}_{000}+{R}_{001}}{2},\\ &&{L}_{000}^{\prime }{X}_{00{ss}^{\prime} }^{\prime }{R}_{{ss}^{\prime} u}^{\prime }={L}_{000}\displaystyle \sum _{t^{\prime} ,t}{P}_{(s^{\prime} ,u),(t^{\prime} ,t)}^{{\rm{R}}}{X}_{00\chi (0{st}^{\prime} )t^{\prime} }{R}_{\chi (0{st}^{\prime} )t^{\prime} t},\\ &&{L}_{v^{\prime} {ss}^{\prime} }{X}_{{ss}^{\prime} 00}{R}_{000}=\displaystyle \sum _{t^{\prime} ,t}{P}_{(v^{\prime} ,s),(t^{\prime} ,t)}^{{\rm{L}}}{L}_{t^{\prime} t\chi ({ts}^{\prime} 0)}^{\prime }{X}_{t\chi ({ts}^{\prime} 0)00}^{\prime }{R}_{000}^{\prime },\\ &&{L}_{000}{X}_{0000}{X}_{00{ss}^{\prime} }{R}_{{ss}^{\prime} u}=\displaystyle \frac{{L}_{000}^{\prime }+{L}_{100}^{\prime }}{2}{X}_{000\chi (00s)}^{\prime }{X}_{0\chi (00s)s\chi ({ss}^{\prime} u)}^{\prime }{R}_{s\chi ({ss}^{\prime} u)u}^{\prime }.\end{eqnarray} \tag{ 18 }$$

The total number of $2\times 16+4\times 8-4=60$ unknowns can be further reduced by exploring the following gauge symmetry

$$\begin{eqnarray}{X}_{{ss}^{\prime} {tt}^{\prime} } & \longrightarrow & {f}_{{ss}^{\prime} }{X}_{{ss}^{\prime} {tt}^{\prime} }{f}_{{tt}^{\prime} }^{-1},\\ {L}_{s^{\prime} {tt}^{\prime} } & \longrightarrow & {L}_{s^{\prime} {tt}^{\prime} }{f}_{{tt}^{\prime} }^{-1},\\ {R}_{{ss}^{\prime} t} & \longrightarrow & {f}_{{ss}^{\prime} }{R}_{{ss}^{\prime} t},\\ {X}_{{ss}^{\prime} {tt}^{\prime} }^{\prime } & \longrightarrow & {g}_{{ss}^{\prime} }{X}_{{ss}^{\prime} {tt}^{\prime} }^{\prime }{g}_{{tt}^{\prime} }^{-1},\\ {L}_{s^{\prime} {tt}^{\prime} }^{\prime } & \longrightarrow & {L}_{s^{\prime} {tt}^{\prime} }^{\prime }{g}_{{tt}^{\prime} }^{-1},\\ {R}_{{ss}^{\prime} t}^{\prime } & \longrightarrow & {g}_{{ss}^{\prime} }{R}_{{ss}^{\prime} t}^{\prime },\end{eqnarray} \tag{ 19 }$$

which conserves the patch ansatz (16), as well as the defining equations (18) for arbitrary nonzero gauge 'fields' ${f}_{{ss}^{\prime} },{g}_{{ss}^{\prime} }$. We can uniquely fix ${f}_{{ss}^{\prime} },{g}_{{ss}^{\prime} }$ by choosing the following gauge ${X}_{00{ss}^{\prime} }={X}_{00{ss}^{\prime} }^{\prime }=1,\forall s,s^{\prime} \in \{0,1\}.$ Numerical experiments suggest some further symmetries which finally inspire the following ansatz

$$\begin{eqnarray}{X}_{0000} & = & 1,\quad & {X}_{0000}^{\prime }=1,\qquad \quad & {L}_{000}=1, & \quad {L}_{000}^{\prime }=1,\\ {X}_{0001} & = & 1,\quad & {X}_{0001}^{\prime }=1,\qquad \quad & {L}_{001}=1, & \quad {L}_{001}^{\prime }={x}_{6},\\ {X}_{0010} & = & 1,\quad & {X}_{0010}^{\prime }=1,\qquad \quad & {L}_{001}=1, & \quad {L}_{001}^{\prime }={x}_{6},\\ {X}_{0011} & = & 1,\quad & {X}_{0011}^{\prime }=1,\qquad \quad & {L}_{010}={x}_{7}, & \quad {L}_{010}^{\prime }={x}_{7},\\ {X}_{0100} & = & {x}_{1},\quad & {X}_{0100}^{\prime }={x}_{1},\qquad \quad & {L}_{011}={x}_{6}, & \quad {L}_{011}^{\prime }=1,\\ {X}_{0101} & = & {x}_{1},\quad & {X}_{0101}^{\prime }={x}_{1},\qquad \quad & {L}_{100}=1, & \quad {L}_{100}^{\prime }={x}_{8},\\ {X}_{0110} & = & {x}_{2},\quad & {X}_{0110}^{\prime }={x}_{3},\qquad \quad & {L}_{101}=1, & \quad {L}_{101}^{\prime }={x}_{9},\\ {X}_{0111} & = & {x}_{4},\quad & {X}_{0111}^{\prime }={x}_{5},\qquad \quad & {L}_{111}={x}_{9}, & \quad {L}_{111}^{\prime }=1,\\ {X}_{1000} & = & {x}_{1},\quad & {X}_{1000}^{\prime }={x}_{1},\qquad \quad & {R}_{000}={x}_{12}, & \quad {R}_{000}^{\prime }=1,\\ {X}_{1001} & = & {x}_{1},\quad & {X}_{1001}^{\prime }={x}_{1},\qquad \quad & {R}_{001}={x}_{13}, & \quad {R}_{001}^{\prime }=1,\\ {X}_{1010} & = & {x}_{1},\quad & {X}_{1010}^{\prime }={x}_{1},\qquad \quad & {R}_{010}={x}_{14}, & \quad {R}_{010}^{\prime }={x}_{15},\\ {X}_{1011} & = & {x}_{1},\quad & {X}_{1011}^{\prime }={x}_{1},\qquad \quad & {R}_{011}={x}_{12}, & \quad {R}_{011}^{\prime }={x}_{16},\\ {X}_{1100} & = & {x}_{5},\quad & {X}_{1100}^{\prime }={x}_{4},\qquad \quad & {R}_{100}={x}_{17}, & \quad {R}_{100}^{\prime }={x}_{18},\\ {X}_{1101} & = & {x}_{5},\quad & {X}_{1101}^{\prime }={x}_{4},\qquad \quad & {R}_{101}={x}_{19}, & \quad {R}_{101}^{\prime }={x}_{18},\\ {X}_{1110} & = & 1,\quad & {X}_{1110}^{\prime }=1,\qquad \quad & {R}_{110}={x}_{20}, & \quad {R}_{110}^{\prime }={x}_{21},\\ {X}_{1111} & = & {x}_{1},\quad & {X}_{1111}^{\prime }={x}_{1},\qquad \quad & {R}_{111}={x}_{20}, & \quad {R}_{111}^{\prime }={x}_{22},\end{eqnarray} \tag{ 20 }$$

with 22 unknown parameters/variables $\{{x}_{i};i=1,\ldots ,22\}$. Assuming that all components are nonvanishing, i.e., ${x}_{j}\ne 0$, the defining relations (18) are equivalent to the following set of polynomial equations

$$\begin{eqnarray}&&{x}_{1}{x}_{2}-{x}_{4}=0,\qquad {x}_{3}{x}_{4}-1=0,\qquad {x}_{1}{x}_{3}-{x}_{5}=0,\qquad {x}_{4}{x}_{5}-{x}_{1}=0,\\ &&{x}_{8}-2{x}_{12}+1=0,\quad {x}_{6}+{x}_{9}-2{x}_{12}=0,\quad {x}_{12}+{x}_{13}-2=0,\\ &&2{x}_{1}{x}_{11}-{x}_{12}-{x}_{13}=0,\ \ 2{x}_{4}-{x}_{1}{x}_{2}({x}_{12}+{x}_{13})=0,\ \ 2{x}_{8}-{x}_{1}{x}_{10}({x}_{12}+{x}_{13})=0,\\ &&2{x}_{14}-(1+{x}_{8}){x}_{15}=0,\quad 2{x}_{13}-(1+{x}_{8}){x}_{16}=0,\quad {x}_{17}-2{x}_{18}+{x}_{19}=0,\\ &&2{x}_{20}-{x}_{3}(1+{x}_{8}){x}_{18}=0,\ \ 2{x}_{19}-{x}_{5}(1+{x}_{8}){x}_{22}=0,\ \ 2{x}_{17}-{x}_{5}(1+{x}_{8}){x}_{21}=0,\\ &&{x}_{1}(\alpha {x}_{7}+(1-\beta ){x}_{11})-{x}_{5}{x}_{6}{x}_{12}=0,\qquad (\beta -\alpha -1){x}_{4}+{x}_{1}{x}_{7}{x}_{12}=0,\\ &&{x}_{1}((1-\alpha ){x}_{7}+\beta {x}_{11})-{x}_{5}{x}_{9}{x}_{12}=0,\qquad (\alpha -\beta -1){x}_{4}+{x}_{1}{x}_{10}{x}_{12}=0,\\ &&(1-\delta ){x}_{12}+\gamma {x}_{14}-{x}_{21}=0,\qquad {x}_{16}+(\gamma -\delta -1){x}_{20}=0,\\ &&\delta {x}_{12}+(1-\gamma ){x}_{14}-{x}_{22}=0,\qquad {x}_{15}+(\delta -\gamma -1){x}_{20}=0.\end{eqnarray} \tag{ 21 }$$

Luckily, this system of nonlinear equations admits a simple solution, which can be compactly written introducing the difference driving parameters:

$$\begin{eqnarray}&&\lambda =\alpha -\beta ,\qquad \mu =\gamma -\delta ,\end{eqnarray} \tag{ 22 }$$

namely:

$$\begin{eqnarray*}{x}_{1} & = & \displaystyle \frac{(\lambda +2)(\mu +2)}{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},{x}_{2}=-\displaystyle \frac{{(\lambda \mu +\lambda -2)}^{2}}{(\lambda +2)(\lambda \mu +\mu -2)},\\ {x}_{3} & = & -\displaystyle \frac{{(\lambda \mu +\mu -2)}^{2}}{(\mu +2)(\lambda \mu +\lambda -2)},\\ {x}_{4} & = & -\displaystyle \frac{(\mu +2)(\lambda \mu +\lambda -2)}{{(\lambda \mu +\mu -2)}^{2}},\\ {x}_{5} & = & -\displaystyle \frac{(\lambda +2)(\lambda \mu +\mu -2)}{{(\lambda \mu +\lambda -2)}^{2}},\\ {x}_{6} & = & -\displaystyle \frac{(\lambda \mu +\lambda -2)((\lambda +1)\mu (\lambda +1-2\alpha )-2(\alpha +1)\lambda -2)}{(\lambda +2)(\lambda \mu +\mu -2)},\\ {x}_{7} & = & -\displaystyle \frac{(\lambda +1)(\lambda \mu +\lambda -2)}{\lambda +2},\\ {x}_{8} & = & \displaystyle \frac{(\lambda -1)(\mu +2)}{\lambda \mu +\mu -2},\\ {x}_{9} & = & \displaystyle \frac{(\lambda \mu +\lambda -2)((\lambda +1)\mu (\lambda +1-2\alpha )-2\alpha \lambda +2)}{(\lambda +2)(\lambda \mu +\mu -2)},\\ {x}_{10} & = & \displaystyle \frac{(\lambda -1)(\lambda \mu +\lambda -2)}{\lambda +2},\end{eqnarray*}$$

$$\begin{eqnarray*}{x}_{11} & = & \displaystyle \frac{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)}{(\lambda +2)(\mu +2)},\\ {x}_{12} & = & \displaystyle \frac{\lambda \mu +\lambda -2}{\lambda \mu +\mu -2},\\ {x}_{13} & = & \displaystyle \frac{(\lambda +2)(\mu -1)}{\lambda \mu +\mu -2},\\ {x}_{14} & = & -\displaystyle \frac{(\lambda +2)(\mu +1)}{\lambda \mu +\mu -2},\\ {x}_{15} & = & -\displaystyle \frac{(\lambda +2)(\mu +1)}{\lambda \mu +\lambda -2},\\ {x}_{16} & = & \displaystyle \frac{(\lambda +2)(\mu -1)}{\lambda \mu +\lambda -2},\\ {x}_{17} & = & -\displaystyle \frac{(\lambda +2)(\lambda (\mu +1)(\mu +1-2\gamma )-2(\gamma +1)\mu -2)}{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},\\ {x}_{18} & = & \displaystyle \frac{(\lambda +2)(\mu +2)}{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},\\ {x}_{19} & = & \displaystyle \frac{(\lambda +2)(\lambda (\mu +1)(\mu +1-2\gamma )-2\gamma \mu +2)}{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},\\ {x}_{20} & = & -\displaystyle \frac{\lambda +2}{\lambda \mu +\lambda -2},\end{eqnarray*}$$

$$\begin{eqnarray}{x}_{21} & = & \displaystyle \frac{\lambda (\mu +1)(\mu +1-2\gamma )-2(\gamma +1)\mu -2}{\lambda \mu +\mu -2},\\ {x}_{22} & = & -\displaystyle \frac{\lambda (\mu +1)(\mu +1-2\gamma )-2\gamma \mu +2}{\lambda \mu +\mu -2}.\end{eqnarray} \tag{ 23 }$$

Note that tensors X and $X^{\prime}$ (components ${x}_{1},\ldots ,{x}_{5}$), which determine the bulk properties of NESS, depend only on the difference parameters $\lambda ,\mu$, while some components of the boundary tensors $L,R,L^{\prime} ,R^{\prime}$ (namely ${x}_{6},{x}_{9},{x}_{17},{x}_{19},{x}_{21},{x}_{22}$) depend explicitly also on the offset parameters $\alpha ,\gamma$. Furthermore, all components are strictly positive, ${x}_{j}\gt 0,j=1,\ldots ,22$, on the open physical domain $(\alpha ,\beta ,\gamma ,\delta )\in {(\mathrm{0,1})}^{\times 4}$. Notice as well that the Gröbner basis algorithm implemented within Mathematica yielded one more solution to the system of equations (21), which however can be excluded by verifying components of the fixed point condition (11) beyond the three-cluster and four-cluster configurations.

(ii) Having explicit expressions for the patch tensors $X,X^{\prime} ,L,L^{\prime} ,R,R^{\prime}$ (23) one can now check that the ansatz (16) solves the full set of component-wise equations of the fixed point condition (11) for a small system size, say for n = 8 where the PSA contains products of two subsequent X ($X^{\prime}$) tensors. This is readily confirmed using computer algebra. Similarly, one may verify the following set of remarkable composition identities

$$\begin{eqnarray}&&\displaystyle \frac{{X}_{{ss}^{\prime} {tt}^{\prime} }^{\prime }{X}_{{tt}^{\prime} {uu}^{\prime} }^{\prime }}{{X}_{{ss}^{\prime} {uu}^{\prime} }^{\prime }}=\displaystyle \frac{{X}_{\chi (v^{\prime} {ss}^{\prime} )s^{\prime} \chi (s^{\prime} {tt}^{\prime} )t^{\prime} }{X}_{\chi (s^{\prime} {tt}^{\prime} )t^{\prime} \chi (t^{\prime} {uu}^{\prime} )u^{\prime} }{X}_{\chi (t^{\prime} {uu}^{\prime} )u^{\prime} \chi (u^{\prime} {zz}^{\prime} )z^{\prime} }}{{X}_{\chi (v^{\prime} {ss}^{\prime} )s^{\prime} \chi (s^{\prime} {uu}^{\prime} )u^{\prime} }{X}_{\chi (s^{\prime} {uu}^{\prime} )u^{\prime} \chi (u^{\prime} {zz}^{\prime} )z^{\prime} }},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{X}_{{ss}^{\prime} {tt}^{\prime} }{X}_{{tt}^{\prime} {uu}^{\prime} }}{{X}_{{ss}^{\prime} {uu}^{\prime} }}=\displaystyle \frac{{X}_{z\chi ({zz}^{\prime} s)s\chi ({ss}^{\prime} t)}^{\prime }{X}_{s\chi ({ss}^{\prime} t)t\chi ({tt}^{\prime} u)}^{\prime }{X}_{t\chi ({tt}^{\prime} u)u\chi ({uu}^{\prime} v)}^{\prime }}{{X}_{z\chi ({zz}^{\prime} s)s\chi ({ss}^{\prime} u)}^{\prime }{X}_{s\chi ({ss}^{\prime} u)u\chi ({uu}^{\prime} v)}^{\prime }},\end{eqnarray} \tag{ 25 }$$

for an arbitrary configuration of indices $s,s^{\prime} ,t,t^{\prime} ,u,u^{\prime} ,z,z^{\prime}$, and $v^{\prime}$ for (24) or v for (25), that means in total $2\times {2}^{9}=1024$ identities.

Now we can make an inductive step. We assume that (16) solves (11) for some even n, writing the PSA compactly as

$$\begin{eqnarray}&&{p}_{\ldots v^{\prime} {ss}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }\;=\;{X}_{\ldots v^{\prime} {ss}^{\prime} }^{{\rm{L}}}{X}_{{ss}^{\prime} {uu}^{\prime} }{X}_{{uu}^{\prime} {zz}^{\prime} }{X}_{{zz}^{\prime} \ldots }^{{\rm{R}}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{p}_{\ldots {zz}^{\prime} {ss}^{\prime} {uu}^{\prime} v\ldots }^{\prime }\;=\;{X}_{\ldots {zz}^{\prime} }^{\prime {\rm{L}}}{X}_{{zz}^{\prime} {ss}^{\prime} }^{\prime }{X}_{{ss}^{\prime} {uu}^{\prime} }^{\prime }{X}_{{uu}^{\prime} v\ldots }^{\prime {\rm{R}}},\end{eqnarray} \tag{ 27 }$$

where ${X}_{\underline{s}}^{{\rm{L}}}$, or ${X}_{\underline{s}}^{\prime {\rm{L}}}$, denote a patch product⁶ of a tensor L, or $L^{\prime}$, with an appropriate number (which could also be zero) of tensors X, or $X^{\prime}$, while ${X}_{\underline{s}}^{{\rm{R}}}$, or ${X}_{\underline{s}}^{\prime {\rm{R}}}$, denote a patch product of an appropriate number of tensors X, or $X^{\prime}$, with R, or $R^{\prime}$. Then, the condition that the PSA

$$\begin{eqnarray}&&{p}_{\ldots v^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }\ =\ {X}_{\ldots v^{\prime} {ss}^{\prime} }^{{\rm{L}}}{X}_{{ss}^{\prime} {tt}^{\prime} }{X}_{{tt}^{\prime} {uu}^{\prime} }{X}_{{uu}^{\prime} {zz}^{\prime} }{X}_{{zz}^{\prime} \ldots }^{{\rm{R}}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{p}_{\ldots {zz}^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} v\ldots }^{\prime }={X}_{\ldots {zz}^{\prime} }^{\prime {\rm{L}}}{X}_{{zz}^{\prime} {ss}^{\prime} }^{\prime }{X}_{{ss}^{\prime} {tt}^{\prime} }^{\prime }{X}_{{tt}^{\prime} {uu}^{\prime} }{X}_{{uu}^{\prime} v\ldots }^{\prime {\rm{R}}},\end{eqnarray} \tag{ 29 }$$

for the system of size $n+2$ again solves the fixed point equations (11), namely

$$\begin{eqnarray}&&\displaystyle \frac{{({U}_{{\rm{e}}}\underline{p})}_{\ldots v^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }}{{({U}_{{\rm{e}}}\underline{p})}_{\ldots v^{\prime} {ss}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }}=\displaystyle \frac{{p}_{\ldots v^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }^{\prime }}{{p}_{\ldots v^{\prime} {ss}^{\prime} {uu}^{\prime} {zz}^{\prime} \ldots }^{\prime }},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{({U}_{{\rm{o}}}\underline{p}^{\prime} )}_{\ldots {zz}^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} v\ldots }}{{({U}_{{\rm{o}}}\underline{p}^{\prime} )}_{\ldots {zz}^{\prime} {ss}^{\prime} {uu}^{\prime} v\ldots }}=\displaystyle \frac{{p}_{\ldots {zz}^{\prime} {ss}^{\prime} {tt}^{\prime} {uu}^{\prime} v\ldots }}{{p}_{\ldots {zz}^{\prime} {ss}^{\prime} {uu}^{\prime} v\ldots }},\end{eqnarray} \tag{ 31 }$$

amounts exactly to X-system (24) and (25).□

We could have even started the induction at n = 6 as an additional set of $2\times {2}^{7}$ composition identities involving the boundary tensors can be straightforwardly verified using computer algebra as well

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{s^{\prime} {tt}^{\prime} }^{\prime }{X}_{{tt}^{\prime} {uu}^{\prime} }^{\prime }}{{L}_{s^{\prime} {uu}^{\prime} }^{\prime }}=\displaystyle \frac{{L}_{s^{\prime} \chi (s^{\prime} {tt}^{\prime} )t^{\prime} }{X}_{\chi (s^{\prime} {tt}^{\prime} )t^{\prime} \chi (t^{\prime} {uu}^{\prime} )u^{\prime} }{X}_{\chi (t^{\prime} {uu}^{\prime} )u^{\prime} \chi (u^{\prime} {zz}^{\prime} )z^{\prime} }}{{L}_{s^{\prime} \chi (s^{\prime} {uu}^{\prime} )u^{\prime} }{X}_{\chi (s^{\prime} {uu}^{\prime} )u^{\prime} \chi (u^{\prime} {zz}^{\prime} )z^{\prime} }},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{X}_{{ss}^{\prime} {tt}^{\prime} }{R}_{{tt}^{\prime} u}}{{R}_{{ss}^{\prime} u}}=\displaystyle \frac{{X}_{z\chi ({zz}^{\prime} s)s\chi ({ss}^{\prime} t)}^{\prime }{X}_{s\chi ({ss}^{\prime} t)t\chi ({tt}^{\prime} u)}^{\prime }{R}_{t\chi ({tt}^{\prime} u)u}^{\prime }}{{X}_{z\chi ({zz}^{\prime} s)s\chi ({ss}^{\prime} u)}^{\prime }{R}_{s\chi ({ss}^{\prime} u)u}^{\prime }}.\end{eqnarray} \tag{ 33 }$$

Having an explicit result for the NESS probability vectors $\underline{p},\underline{p}^{\prime}$ (16), (20) and (23) at hand we shall now address a natural question of computation of physical observables, such as density profiles and density–density correlations in the steady state. To start with, let us define the nonequilibrium partition function via normalization of the total probability of the state, and the corresponding transfer matrix.

We observe the following obvious identities, noting again that $n=2m$,

$$\begin{eqnarray}&&{Z}_{n}=\displaystyle \sum _{\underline{s}\in { \mathcal C }}\;{p}_{\underline{s}}=\underline{l}\cdot {T}^{m-2}\underline{r},\qquad {Z}_{n}^{\prime }=\displaystyle \sum _{\underline{s}\in { \mathcal C }}\;{p}_{\underline{s}}^{\prime }=\underline{l}^{\prime} \cdot T{^{\prime} }^{m-2}\underline{r}^{\prime} ,\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&\underline{l}={\underline{l}}_{0}+{\underline{l}}_{1},\quad \underline{l}^{\prime} ={\underline{l}}_{0}^{\prime }+{\underline{l}}_{1}^{\prime },\quad \underline{r}={\underline{r}}_{0}+{\underline{r}}_{1},\quad \underline{r}^{\prime} ={\underline{r}}_{0}^{\prime }+{\underline{r}}_{1}^{\prime },\end{eqnarray} \tag{ 35 }$$

and ${\underline{l}}_{u},{\underline{l}}_{u}^{\prime },{\underline{r}}_{u},{\underline{r}}_{u}^{\prime },u\in \{0,1\}$, are vectors from ${{\mathbb{R}}}^{4}$ with components labelled as $0,1,2,3$, namely:

$$\begin{eqnarray}{({\underline{l}}_{u^{\prime} })}_{2s+s^{\prime} } & = & {L}_{u^{\prime} {ss}^{\prime} },\quad {({\underline{r}}_{u})}_{2s+s^{\prime} }={R}_{{ss}^{\prime} u},\\ {({\underline{l}}_{u^{\prime} }^{\prime })}_{2s+s^{\prime} } & = & {L}_{u^{\prime} {ss}^{\prime} }^{\prime },\quad {({\underline{r}}_{u}^{\prime })}_{2s+s^{\prime} }={R}_{{ss}^{\prime} u}^{\prime },\end{eqnarray} \tag{ 36 }$$

and $T,T^{\prime} \in {\rm{End}}({{\mathbb{R}}}^{4})$ are 4 × 4 transfer matrices with components

$$\begin{eqnarray}&&{T}_{2s+s^{\prime} ,2u+u^{\prime} }={X}_{{ss}^{\prime} {uu}^{\prime} },\quad {T}_{2s+s^{\prime} ,2u+u^{\prime} }^{\prime }={X}_{{ss}^{\prime} {uu}^{\prime} }^{\prime }.\end{eqnarray} \tag{ 37 }$$

Straightforward calculation by means of the explicit solution (23) shows that T and $T^{\prime}$ are similar, i.e. $\exists W\in {\rm{End}}({{\mathbb{R}}}^{4})$,

$$\begin{eqnarray}&&T^{\prime} ={{WTW}}^{-1},\end{eqnarray} \tag{ 38 }$$

such that also

$$\begin{eqnarray}&&\underline{l}^{\prime} =\kappa \underline{l}{W}^{-1},\quad \underline{r}^{\prime} ={\kappa }^{-1}W\underline{r},\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{(\lambda +2)(\lambda \mu +\mu -2)}{(\mu +2)(\lambda \mu +\lambda -2)},\end{eqnarray} \tag{ 40 }$$

and consequently

$$\begin{eqnarray}&&{Z}_{n}={Z}_{n}^{\prime },\end{eqnarray} \tag{ 41 }$$

for any pair of driving parameters $\lambda ,\mu \in (-1,1)$, which is nothing but the conservation of total probability.

Writing a pair of independent spectral parameters as

$$\begin{eqnarray}&&\omega \equiv {x}_{4}=\varphi (\lambda ,\mu ),\ \ \xi \equiv {x}_{5}=\varphi (\mu ,\lambda ),\quad \varphi (\lambda ,\mu )\equiv -\displaystyle \frac{(\mu +2)(\lambda \mu +\lambda -2)}{{(\lambda \mu +\mu -2)}^{2}}\end{eqnarray} \tag{ 42 }$$

we find compact expressions for the transfer- and the intertwining matrices

$$\begin{eqnarray}T=\left(\begin{array}{cccc}1 & 1 & 1 & 1\\ \xi \omega & \xi \omega & \frac{1}{\xi } & \omega \\ \xi \omega & \xi \omega & \xi \omega & \xi \omega \\ \xi & \xi & 1 & \xi \omega \end{array}\right),\quad T^{\prime} =\left(\begin{array}{cccc}1 & 1 & 1 & 1\\ \xi \omega & \xi \omega & \frac{1}{\omega } & \xi \\ \xi \omega & \xi \omega & \xi \omega & \xi \omega \\ \omega & \omega & 1 & \xi \omega \end{array}\right),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}W=\left(\begin{array}{cccc}\frac{\omega +1}{\omega } & \frac{1}{\omega } & -\frac{1}{\xi \omega } & -\frac{1}{\xi \omega }\\ -\frac{\omega +1}{\omega } & -\frac{1}{\omega } & \frac{\xi +1}{\xi \omega } & \frac{\xi +1}{\xi \omega }\\ \xi & \xi & 0 & -1\\ 0 & 0 & 0 & 1\end{array}\right).\end{eqnarray} \tag{ 44 }$$

Remarkably, T and $T^{\prime}$ are swapped upon an exchange of driving parameters λ and μ, or equivalently, exchanging the spectral parameters ω and ξ:

$$\begin{eqnarray}&&T(\omega ,\xi )\equiv T^{\prime} (\xi ,\omega ).\end{eqnarray} \tag{ 45 }$$

In fact, W acts as an intertwiner connecting two subsequent temporal layers

$$\begin{eqnarray}&&W(\omega ,\xi )T(\omega ,\xi )=T(\xi ,\omega )W(\omega ,\xi ),\end{eqnarray} \tag{ 46 }$$

satisfying an inversion identity

$$\begin{eqnarray}&&{[W(\omega ,\xi )]}^{-1}=W(\xi ,\omega ).\end{eqnarray} \tag{ 47 }$$

On the other hand, equation (38) or (46) can be interpreted as an isospectral problem (or discrete space–time zero-curvature condition) with a Lax pair $\{T(\omega ,\xi ),W(\omega ,\xi )\}$. This is a clear indication of Lax integrability of our nonequilibrium steady state.

The transfer matrices $T,T^{\prime}$ are singular with, in general, three nonvanishing eigenvalues, $T=V{\rm{diag}}({\tau }_{1},{\tau }_{2},{\tau }_{3},0){V}^{-1},T^{\prime} =V^{\prime} {\rm{diag}}({\tau }_{1},{\tau }_{2},{\tau }_{3},0)V{^{\prime} }^{-1}$, with $W=V^{\prime} {V}^{-1}$, reading explicitly

$$\begin{eqnarray}{\tau }_{1} & = & \displaystyle \frac{{(\lambda \mu -4)}^{2}}{(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},\\ {\tau }_{\mathrm{2,3}} & = & \displaystyle \frac{\lambda \mu (\lambda +\mu +8)+4(\lambda +\mu )\pm (\lambda -\mu )\sqrt{D}}{2(\lambda \mu +\lambda -2)(\lambda \mu +\mu -2)},\\ D & = & \lambda \mu (\lambda \mu -12)-8(\lambda +\mu ).\end{eqnarray} \tag{ 48 }$$

We remark two important facts: (i) in the whole open parameter domain $\mu ,\lambda \in (-1,1),{\tau }_{1}$ represents the leading eigenvalue ${\tau }_{1}\gt | {\tau }_{\mathrm{2,3}}|$. (ii) In the limit $\mu ,\lambda \to 1$, the three eigenvalues collapse, ${\tau }_{1}/{\tau }_{\mathrm{2,3}}\to 1$, so we should see long range correlations there (to be discussed below); namely, the correlation length should diverge as $\lambda ,\mu \to 1$. Also, an algebraic curve along which the discriminant vanishes $D(\lambda ,\mu )=0$ is potentially interesting, as it signals discontinuous behavior separating the regime of $| {\tau }_{2}| =| {\tau }_{3}|$ from the regime of $| {\tau }_{2}| \ne | {\tau }_{3}|$.

A key feature of integrability of our boundary driven CA is a compatibility condition between the bulk transfer matrix and the boundary vectors, namely

$$\begin{eqnarray}&&\underline{l}T={\tau }_{1}\underline{l},\quad \underline{l}^{\prime} T^{\prime} ={\tau }_{1}\underline{l}^{\prime} ,\quad T\underline{r}={\tau }_{1}\underline{r},\quad T^{\prime} \underline{r}^{\prime} ={\tau }_{1}\underline{r}^{\prime} ,\end{eqnarray} \tag{ 49 }$$

relations, which can be readily verified from our analytic solution. This yields a particularly simple expression for the partition sum

$$\begin{eqnarray}&&{Z}_{n}=(\underline{l}\cdot \underline{r}){\tau }_{1}^{m-2}.\end{eqnarray} \tag{ 50 }$$

We are now ready to compute the density profiles. Let us define the following pair of diagonal matrices

$$\begin{eqnarray}&&{D}_{{\rm{e}}}={\rm{diag}}(0,0,1,1),\quad {D}_{{\rm{o}}}={\rm{diag}}(0,1,0,1).\end{eqnarray} \tag{ 51 }$$

Then the steady-state density profiles on even and odd bulk sites express as:

$$\begin{eqnarray}{\rho }_{2k} & = & \displaystyle \frac{1}{{Z}_{n}}\displaystyle \sum _{{s}_{1},{s}_{2},\ldots {s}_{n}}{s}_{2k}\;{p}_{{s}_{1},{s}_{2},\ldots {s}_{n}}=\displaystyle \frac{1}{{Z}_{n}}\underline{l}\cdot {T}^{k-1}{D}_{{\rm{e}}}{T}^{m-1-k}\underline{r},\\ {\rho }_{2k+1} & = & \displaystyle \frac{1}{{Z}_{n}}\displaystyle \sum _{{s}_{1},{s}_{2},\ldots {s}_{n}}{s}_{2k+1}\;{p}_{{s}_{1},{s}_{2},\ldots {s}_{n}}=\displaystyle \frac{1}{{Z}_{n}}\underline{l}\cdot {T}^{k-1}{D}_{{\rm{o}}}{T}^{m-1-k}\underline{r},\end{eqnarray} \tag{ 52 }$$

for $k=1,2,\ldots \;m-1$, while at the boundary sites

$$\begin{eqnarray}&&{\rho }_{1}=\displaystyle \frac{1}{{Z}_{n}}{\underline{l}}_{1}\cdot {T}^{m-2}\underline{r},\quad {\rho }_{n}=\displaystyle \frac{1}{{Z}_{n}}\underline{l}\cdot {T}^{m-2}{\underline{r}}_{1}.\end{eqnarray} \tag{ 53 }$$

The compatibility conditions (49) immediately imply flat—ballistic—density profiles apart from the boundary sites

$$\begin{eqnarray}{\rho }_{1} & = & \displaystyle \frac{{\underline{l}}_{1}\cdot \underline{r}}{\underline{l}\cdot \underline{r}}=\displaystyle \frac{2(1-\alpha )(\lambda \mu +\lambda +\mu )+(\lambda -2)(\lambda \mu -4)}{2(\lambda +\mu +8-\lambda \mu )},\\ {\rho }_{j} & = & \displaystyle \frac{\underline{l}\cdot {D}_{{\rm{e}}}\underline{r}}{\underline{l}\cdot \underline{r}}=\displaystyle \frac{\underline{l}\cdot {D}_{{\rm{o}}}\underline{r}}{\underline{l}\cdot \underline{r}}=\displaystyle \frac{\lambda +\mu +4}{\lambda +\mu +8-\lambda \mu },\qquad {\rm{for}}\ \ 1\lt j\lt n,\\ {\rho }_{n} & = & \displaystyle \frac{\underline{l}\cdot {\underline{r}}_{1}}{\underline{l}\cdot \underline{r}}=\displaystyle \frac{2(1-\gamma )(\lambda \mu +\lambda +\mu )+(\mu -2)(\lambda \mu -4)}{2(\lambda +\mu +8-\lambda \mu )}.\end{eqnarray} \tag{ 54 }$$

Interestingly, the bulk steady-state density can only take values in the interval ${\rho }_{j}\in \left(\frac{2}{5},\frac{2}{3}\right)$, with the extreme value of maximum density $\frac{2}{3}$ (minimum density $\frac{2}{5}$) reached in the limit $\mu ,\lambda \to 1$ ($\mu ,\lambda \to -1$).

Similarly one computes a steady-state two-point density–density correlation function, say for even–even sites $2k,2k^{\prime}$, assuming $k\lt k^{\prime}$:

$$\begin{eqnarray}{\tilde{C}}_{2k,2k^{\prime} } & = & \displaystyle \frac{1}{{Z}_{n}}\displaystyle \sum _{{s}_{1},{s}_{2},\ldots {s}_{n}}{s}_{2k}{s}_{2k^{\prime} }\;{p}_{{s}_{1},{s}_{2},\ldots {s}_{n}}=\displaystyle \frac{1}{{Z}_{n}}\underline{l}\cdot {T}^{k-1}{D}_{{\rm{e}}}{T}^{k^{\prime} -k}{D}_{{\rm{e}}}{T}^{m-1-k^{\prime} }\underline{r}\\ & = & \displaystyle \frac{{D}_{{\rm{e}}}\underline{l}\cdot {T}^{k^{\prime} -k}{D}_{{\rm{e}}}\underline{r}}{\underline{l}\cdot {T}^{k^{\prime} -k}\underline{r}},\end{eqnarray} \tag{ 55 }$$

and similarly for other pairs of sites, even–odd, odd–even and odd–odd, exchanging the corresponding ${D}_{{\rm{e}}}$ with ${D}_{{\rm{o}}}$. Then, the connected two-point function is defined as

$$\begin{eqnarray}&&{C}_{j,j^{\prime} }={\tilde{C}}_{j,j^{\prime} }-{\rho }_{j}{\rho }_{j^{\prime} }.\end{eqnarray} \tag{ 56 }$$

Writing an eigenvalue decomposition of the transfer matrix as

$$\begin{eqnarray}&&T=\displaystyle \sum _{\nu =1}^{3}{\tau }_{\nu }\;{\underline{\psi }}_{\nu }\otimes {\underline{\phi }}_{\nu }\end{eqnarray} \tag{ 57 }$$

where ${\underline{\phi }}_{1}=\underline{l}/\sqrt{\underline{l}\cdot \underline{r}},{\underline{\psi }}_{1}=\underline{r}/\sqrt{\underline{l}\cdot \underline{r}}$, so that ${\phi }_{\nu }\cdot {\psi }_{\nu ^{\prime} }={\delta }_{\nu ,\nu ^{\prime} }$, we can rewrite the connected correlation explicitly as

$$\begin{eqnarray}&&{C}_{2k,2k^{\prime} }=\displaystyle \frac{({D}_{{\rm{e}}}\underline{l}\cdot {\underline{\psi }}_{2})({\underline{\phi }}_{2}\cdot {D}_{{\rm{e}}}\underline{r})}{\underline{l}\cdot \underline{r}}{\left(\displaystyle \frac{{\tau }_{2}}{{\tau }_{1}}\right)}^{k^{\prime} -k}+\displaystyle \frac{({D}_{{\rm{e}}}\underline{l}\cdot {\underline{\psi }}_{3})({\underline{\phi }}_{3}\cdot {D}_{{\rm{e}}}\underline{r})}{\underline{l}\cdot \underline{r}}{\left(\displaystyle \frac{{\tau }_{3}}{{\tau }_{1}}\right)}^{k^{\prime} -k}.\quad \end{eqnarray} \tag{ 58 }$$

This demonstrates that the connected two-point correlation function (56) in the bulk $1\lt j,j^{\prime} \lt n$ is indeed only a function of the difference of indices (positions) which decays exponentially $\sim \mathrm{exp}(-| j-j^{\prime} | /{\ell })$ with the correlation scale ${\ell }=1/\mathrm{log}| {\tau }_{1}/{\tau }_{2}|$, and depends only on the difference driving parameters $\lambda ,\mu$ and not on $\alpha ,\beta ,\gamma ,\delta$ separately. See figure 5 for an example. Extension of such calculations to higher k-point connected correlations is straightforward; they would all decay exponentially $\sim \mathrm{exp}(-| j-j^{\prime} | /{\ell })$ in difference of any pair of adjacent spatial coordinates $j,j^{\prime}$.

Zoom In Zoom Out Reset image size

And finally, let us compute the steady state soliton currents. The current can be computed as the density of right-movers minus the density of left movers [1]. The density of right-movers, computed as

$$\begin{eqnarray}&&{J}_{{\rm{R}}}={Z}_{n}^{-1}\displaystyle \sum _{\underline{s}}{s}_{2k}{s}_{2k+1}{p}_{\underline{s}}={Z}_{n}^{-1}\underline{l}\cdot {T}^{k-1}{D}_{{\rm{e}}}{D}_{{\rm{o}}}{T}^{m-1-k}\underline{r}=\displaystyle \frac{\underline{l}\cdot {D}_{{\rm{e}}}{D}_{{\rm{o}}}\underline{r}}{\underline{l}\cdot \underline{r}},\end{eqnarray} \tag{ 59 }$$

is independent of location k in the steady state and reads explicitly

$$\begin{eqnarray}&&{J}_{{\rm{R}}}=\displaystyle \frac{\lambda +2}{\lambda +\mu +8-\lambda \mu }.\end{eqnarray} \tag{ 60 }$$

Similarly, the density of left-movers reads

$$\begin{eqnarray}{J}_{{\rm{L}}} & = & {Z}_{n}^{-1}\displaystyle \sum _{\underline{s}}\;{s}_{2k+1}{s}_{2k+2}{p}_{\underline{s}}=\underline{l}\cdot {T}^{k-1}{D}_{{\rm{o}}}{{TD}}_{{\rm{e}}}{T}^{m-1-k}\underline{r}=\displaystyle \frac{\underline{l}\cdot {D}_{{\rm{o}}}{{TD}}_{{\rm{e}}}\underline{r}}{{\tau }_{1}\underline{l}\cdot \underline{r}}\\ & = & \displaystyle \frac{\mu +2}{\lambda +\mu +8-\lambda \mu },\end{eqnarray} \tag{ 61 }$$

with the overall steady state soliton current

$$\begin{eqnarray}&&J={J}_{{\rm{R}}}-{J}_{{\rm{L}}}=\displaystyle \frac{\lambda -\mu }{\lambda +\mu +8-\lambda \mu }.\end{eqnarray} \tag{ 62 }$$

Note the expected linear-response behavior for small driving parameters, namely the current becomes linearly proportional to the bias $J\sim \frac{1}{8}(\lambda -\mu )$.

Let us now try to approach the problem of finding possibly a complete set of independent conservation laws of RCA54 dynamics from a more formal nonequilibrium point of view. We shall take an approach analogous to the construction of quasilocal conservation laws of integrable quantum spin chains via dissipative boundary driving [8, 9]. First, let us equip the space ${ \mathcal S }$ of probability distributions with a Hilbert inner product

$$\begin{eqnarray}&&(\underline{u}| \underline{v})={2}^{-n}\;\underline{u}\cdot \underline{v}={2}^{-n}\displaystyle \sum _{\underline{s}\in { \mathcal C }}{u}_{\underline{s}}{v}_{\underline{s}},\quad \underline{u},\underline{v}\in { \mathcal S }.\end{eqnarray} \tag{ 63 }$$

Writing orthogonal basis vectors on ${{\mathbb{R}}}^{2}$ as ${\underline{\omega }}^{0}=(1,1)$ and ${\underline{\omega }}^{1}=(1,-1)$, we can write a convenient orthonormal basis $\{{\underline{\omega }}^{{b}_{1},{b}_{2},\ldots ,{b}_{n}},{b}_{j}\in {{\mathbb{Z}}}_{2}\}$ of ${ \mathcal S }$ as

$$\begin{eqnarray}&&{\underline{\omega }}^{\underline{b}}={\underline{\omega }}^{{b}_{1}}\otimes {\underline{\omega }}^{{b}_{2}}\otimes \cdots \;\otimes {\underline{\omega }}^{{b}_{n}},\qquad ({\underline{\omega }}^{\underline{b}}| {\underline{\omega }}^{\underline{b}^{\prime} })={\delta }_{\underline{b},\underline{b}^{\prime} },\end{eqnarray} \tag{ 64 }$$

meaning that

$$\begin{eqnarray}&&{\underline{\omega }}_{\underline{s}}^{\underline{b}}=\displaystyle \prod _{j=1}^{n}{\omega }_{{s}_{j}}^{{b}_{j}}={(-1)}^{{\displaystyle \sum }_{j=1}^{n}{b}_{j}{s}_{j}}.\end{eqnarray} \tag{ 65 }$$

The PSA (16) for the un-normalized NESS then immediately translates to a more standard form of a matrix product ansatz

$$\begin{eqnarray}&&\underline{p}=\displaystyle \sum _{\underline{b}\in {\{\mathrm{0,1}\}}^{n}}({\underline{l}}^{{b}_{1}}\cdot {\sigma }_{1}^{{b}_{2}}{\sigma }_{2}^{{b}_{3}}T{\sigma }_{1}^{{b}_{4}}{\sigma }_{2}^{{b}_{5}}T\cdots T{\sigma }_{1}^{{b}_{n-2}}{\sigma }_{2}^{{b}_{n-1}}{\underline{r}}^{{b}_{n}}){\underline{\omega }}^{\underline{b}},\end{eqnarray} \tag{ 66 }$$

and analogous expression for $\underline{p}^{\prime}$, with contravariant boundary vectors

$$\begin{eqnarray}&&{\underline{l}}^{b}=\displaystyle \frac{1}{2}\displaystyle \sum _{s\in \{0,1\}}{(-1)}^{{bs}}{\underline{r}}_{s},\quad {\underline{r}}^{b}=\displaystyle \frac{1}{2}\displaystyle \sum _{s\in \{0,1\}}{(-1)}^{{bs}}{\underline{r}}_{s},\end{eqnarray} \tag{ 67 }$$

and ${\sigma }_{1,2}^{b},b\in \{0,1\}$ are diagonal 4 × 4 matrices, defined as

$$\begin{eqnarray}{\sigma }_{1}^{b}={\sigma }^{b}\otimes {{\mathbb{1}}}_{2},\quad {\sigma }_{2}^{b}={{\mathbb{1}}}_{2}\otimes {\sigma }^{b},\quad {\sigma }^{b}=\displaystyle \frac{1}{2}\left(\begin{array}{cc}1 & 0\\ 0 & {(-1)}^{b}\end{array}\right).\end{eqnarray} \tag{ 68 }$$

Clearly, it can be directly verified that with our definitions (65), (67) and (68), the expression (66) is equivalent to the first line of equations (16), since

$$\begin{eqnarray}\displaystyle \sum _{b}\;{\underline{l}}^{b}{\omega }_{s}^{b}={\underline{l}}_{s},\quad \displaystyle \sum _{b}\;{\underline{r}}^{b}{\omega }_{s}^{b}={\underline{r}}_{s},\quad \displaystyle \sum _{b}\;{\sigma }^{b}{\omega }_{s}^{b}=\left(\begin{array}{cc}{\delta }_{s,0} & 0\\ 0 & {\delta }_{s,1}\end{array}\right).\end{eqnarray} \tag{ 69 }$$

A special state ${({\underline{\omega }}^{0})}^{\otimes n}={\underline{\omega }}^{00\ldots 0}$ represents a uniform distribution over all ${2}^{n}$ configurations ('infinite temperature state'). Any state of the form ${\underline{\psi }}^{(k,r)}={({\underline{\omega }}^{0})}^{\otimes k}\otimes \underline{v}\otimes {({\underline{\omega }}^{0})}^{\otimes (n-k-r)}$ for some ${2}^{r}$ dimensional vectors $\underline{v}\in {({{\mathbb{R}}}^{2})}^{\otimes r}$ is defined as a r-local state, supported on sites $[k,k+r-1]$, and its components may depend only on the coordinates from the supported set ${s}_{k},\ldots ,{s}_{k+r-1}$, namely ${\psi }_{{s}_{1},{s}_{2},\ldots ,{s}_{n}}^{(k,r)}={v}_{{s}_{k},{s}_{k+1},\ldots ,{s}_{k+r-1}}$. In fact, introduction of the Hilbert space metric (63) identifies the state space with its dual—the space of observables, so one may interpret a vector ${\underline{\psi }}^{(k,r)}$ also as a r-local observable. For example, ${\underline{\rho }}^{(j)}={({\underline{\omega }}^{0})}^{\otimes (j-1)}\otimes (0,1)\otimes {({\underline{\omega }}^{0})}^{\otimes (n-j)}$ is the density, with expectation (52) given as

$$\begin{eqnarray}&&{\rho }_{j}=({\underline{\rho }}^{(j)}| \underline{p}).\end{eqnarray} \tag{ 70 }$$

Definition 1. A local conservation law $\underline{Q}\in { \mathcal S }$ of a boundary driven RCA is defined as an extensive sum of a shifted r-local observable (for some even integer r independent of size n), written in terms of a vector $\underline{q}\in {{\mathbb{R}}}^{{2}^{r}}$,

$$\begin{eqnarray}&&\underline{Q}=\displaystyle \sum _{k=1}^{(n-r)/2}{({\underline{\omega }}^{0})}^{\otimes (2k-1)}\otimes \underline{q}\otimes {({\underline{\omega }}^{0})}^{\otimes (n-r-2k+1)}\end{eqnarray} \tag{ 71 }$$

for which the time-difference in one step is localized near the boundaries of the system. More precisely,

$$\begin{eqnarray}&&U\underline{Q}-\underline{Q}=\underline{g}\otimes {({\underline{\omega }}^{0})}^{\otimes (n-r^{\prime} )}+{({\underline{\omega }}^{0})}^{\otimes (n-r^{\prime} )}\otimes \underline{h},\end{eqnarray} \tag{ 72 }$$

for some remainder observables, specified by vectors $\underline{g},\underline{h}\in {{\mathbb{R}}}^{{2}^{r^{\prime} }}$, localized near boundaries with n-independent support size $r^{\prime}$.

Since, when approaching the thermodynamic limit $n\to \infty$, the square norm of translationally invariant sum of local observables is extensive in n,

$$\begin{eqnarray}&&(\underline{Q}| \underline{Q})=\left(\displaystyle \frac{\underline{q}\cdot \underline{q}}{{2}^{r}}\right)\left(\displaystyle \frac{n-r}{2}\right)\propto n,\end{eqnarray} \tag{ 73 }$$

while the remainder—rhs of (72) has a bounded (in n) norm, we can conclude that such $\underline{Q}$ is exactly conserved in the bulk in the thermodynamic limit. A formal proof, invoking causality of RCA in place of the Lieb–Robinson bound, would be a straightforward extension of an analogous result for quantum chains [6].

We shall now derive exact conservation laws of CA rule 54 using exactly the same strategy as for dissipatively boundary driven quantum chains [8, 9]. The NESS probability vector $\underline{p}$ is already a potential candidate for a conservation law, equation (72), since

$$\begin{eqnarray}&&U\underline{p}-\underline{p}=0,\end{eqnarray} \tag{ 74 }$$

provided it could be identified with a local observable. This is possible in the trivial (equlibrium) case of zero biases $\lambda =\mu =0$, and $\alpha ,\gamma =0$, where NESS is trivial ${\underline{p}}^{0}={({\underline{\omega }}^{0})}^{\otimes n}$. Let us set $\alpha =\gamma =0$ for simplicity in the following, while considering arbitrary values of $\alpha ,\gamma$ would only alter the boundary cells and hence the remainder terms $\underline{g},\underline{h}$ and not the bulk density $\underline{q}$. Taking a derivative of equation (74) with expression (66), with respect to, say λ, and putting $\lambda =\mu =0$ at the end, we obtain exactly a conservation law (72), by identifying:

$$\begin{eqnarray}&&U=U{| }_{\lambda =\mu =0},\qquad {\underline{Q}}^{\lambda }={\partial }_{\lambda }\underline{p}{| }_{\lambda =\mu =0},\end{eqnarray} \tag{ 75 }$$

with $(r=4)$-local bulk density

$$\begin{eqnarray}&&{\underline{q}}^{\lambda }=\displaystyle \sum _{a,a^{\prime} ,b,b^{\prime} \in \{0,1\}}(\underline{v}\cdot {\sigma }_{1}^{a}{\sigma }_{2}^{a^{\prime} }{T}_{\lambda }{\sigma }_{1}^{b}{\sigma }_{2}^{b^{\prime} }{T}_{0}\underline{v}){\underline{\omega }}^{{aa}^{\prime} {bb}^{\prime} },\end{eqnarray} \tag{ 76 }$$

where

$$\begin{eqnarray*}\underline{v} & = & {\underline{l}}^{0}{| }_{\lambda =\mu =0}={\underline{r}}^{0}{| }_{\lambda =\mu =0}=(1,1,1,1),\\ {T}_{\lambda } & = & {\partial }_{\lambda }T{| }_{\lambda =\mu =0}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 1 & 1 & -\frac{3}{2} & -\frac{1}{2}\\ 1 & 1 & 1 & 1\\ \frac{3}{2} & \frac{3}{2} & 0 & 1\end{array}\right),\\ {T}_{0} & = & T{| }_{\lambda =\mu =0}=\left(\begin{array}{cccc}1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\end{array}\right).\end{eqnarray*}$$

The local remainder (boundary) terms $\underline{g},\underline{h}$ are generated through terms where ${\partial }_{\lambda }$ either hits ${E}^{\alpha ,\alpha -\lambda }$ in the propagator U, or the boundary vectors ${\underline{l}}^{b}$, or ${\underline{r}}^{b}$ in the expression (66). Similarly, we obtain another local conservation law by differentiating with respect to μ instead of λ:

$$\begin{eqnarray}{\underline{Q}}^{\mu } & = & {\partial }_{\mu }\underline{p}{| }_{\lambda =\mu =0},\\ {\underline{q}}^{\mu } & = & \displaystyle \sum _{a,a^{\prime} ,b,b^{\prime} \in \{0,1\}}4(\underline{v}\cdot {\sigma }_{1}^{a}{\sigma }_{2}^{a^{\prime} }{T}_{\mu }{\sigma }_{1}^{b}{\sigma }_{2}^{b^{\prime} }\underline{v}){\underline{\omega }}^{{aa}^{\prime} {bb}^{\prime} },\\ {T}_{\mu } & = & {\partial }_{\mu }T{| }_{\lambda =\mu =0}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 1 & 1 & \frac{1}{2} & \frac{3}{2}\\ 1 & 1 & 1 & 1\\ -\frac{1}{2} & -\frac{1}{2} & 0 & 1\end{array}\right).\end{eqnarray} \tag{ 77 }$$

Writing ${\underline{Q}}^{\pm }=\frac{1}{2}({\underline{Q}}^{\lambda }\pm {\underline{Q}}^{\mu })$, and consequently ${\underline{q}}^{\pm }=\frac{1}{2}({\underline{q}}^{\lambda }\pm {\underline{q}}^{\mu })$, we find for the density of ${\underline{Q}}^{-}$

$$\begin{eqnarray}&&{\underline{q}}^{-}={\underline{\omega }}^{1100}-{\underline{\omega }}^{0110}+{\underline{\omega }}^{0010}-{\underline{\omega }}^{1000},\end{eqnarray} \tag{ 78 }$$

or in terms of explicit dependence on cell occupation numbers s_j

$$\begin{eqnarray}&&{q}_{{ss}^{\prime} {tt}^{\prime} }^{-}=4({ss}^{\prime} -s^{\prime} t),\end{eqnarray} \tag{ 79 }$$

which is exactly (4 times) the conserved net soliton current (62) as discovered in [1].

The second conservation law is less trivial

$$\begin{eqnarray}&&{q}_{{ss}^{\prime} {tt}^{\prime} }^{+}=4(s+s^{\prime} +{ss}^{\prime} t+s^{\prime} {tt}^{\prime} )-6({ss}^{\prime} +s^{\prime} t),\end{eqnarray} \tag{ 80 }$$

and to best of our knowledge has not been discussed before. We note that both quantities

$$\begin{eqnarray}&&{Q}^{\pm }=\displaystyle \sum _{j}\;{q}_{{s}_{j},{s}_{j+1},{s}_{j+2},{s}_{j+3}}^{\pm }\end{eqnarray} \tag{ 81 }$$

should be exactly conserved for a purely deterministic RCA with periodic boundary conditions.

So far we have studied only the properties of the steady state. An obvious follow up question concerns studying the full relaxation dynamics to the steady state, which amounts to studying the spectrum of decay modes, i.e. eigenvalues $\{{{\rm{\Lambda }}}_{j},j=1,\ldots \;{2}^{n}\}$ of the Markov matrix

$$\begin{eqnarray}&&U{\underline{p}}_{j}={{\rm{\Lambda }}}_{j}{\underline{p}}_{j}.\end{eqnarray} \tag{ 82 }$$

The leading eigenvalue ${{\rm{\Lambda }}}_{1}=1$ corresponds to NESS, while all the others, corresponding to the so-called decay modes, lie strictly inside the unit circle $| {{\rm{\Lambda }}}_{j}| \lt 1,j\gt 1$, as following from our theorem 1. So far we have not been able to provide any exact or rigorous results on the decay modes—and the progress here could be very difficult as we need to devise a particular kind of Bethe ansatz—so we instead report some intriguing results of numerical computations which should strongly motivate further study. Apart from the Markov eigenvalues $\{{{\rm{\Lambda }}}_{j}\}$ we also analyzed the complexity of the corresponding eigenvectors by calculating their Schmidt rank—the number of nonvanishing singular values of the ${2}^{n/2}\times {2}^{n/2}$ matrix ${({P}_{j})}_{\underline{s},\underline{s}^{\prime} }\equiv {({p}_{j})}_{{s}_{1},{s}_{2},\ldots {s}_{n/2},{s}_{1}^{\prime },{s}_{2}^{\prime },\ldots {s}_{n/2}^{\prime }}$. For NESS (j = 1) the Schmidt number is always 3 as we have proven above (e.g. it follows from the fact that the transfer matrices $T,T^{\prime}$ have rank 3, see (48)). Remarkably, there seem to be always two other eigenvectors (decay modes) with Schmidt number 3. The subleading eigenvector (j = 2) corresponds to Schmidt number 6 (independent of n!) which gives hope that this decay mode would be analytically tractable. Curiously, there appears to be even an eigenvector of lower complexity than NESS (Schmidt number 2) with eigenvalue ${\rm{\Lambda }}=-1/2$. The rest of the spectrum is organized in bands with end-points corresponding to nondegenerate eigenvectors with Schmidt number 6. Since a picture says more than a thousand words, the reader is welcome to go and stare at the figure 6.

Zoom In Zoom Out Reset image size

In conclusion, we have presented an exact analytic treatment of the steady-state properties of a strongly interacting deterministic many-body system driven by stochastic boundaries. We considered arguably the simplest possible strongly interacting bulk dynamics which possesses certain features of integrability like solitons, namely the reversible ${{\mathbb{Z}}}_{2}$ cellular automaton with global conservation laws. To facilitate our analysis we have developed a novel algebraic ansatz for describing strongly correlated classical many-body probability states, namely the PSA. We expect that our approach should be applicable for constructing nonequilibrium steady-states of general classical deterministic integrable interacting theories [4] driven by compatible Markovian stochastic boundaries. The fundamental relation proposed here, which needs to be generalized to other integrable models, is a particular 'fusion', or composition formula (24) and (25) which we propose to call the X-system. In an integrable lattice model with a continuous dynamical variable⁷ ${\varphi }_{x}$ at each physical site x, the patch tensor X would be a function of four variables $X({\varphi }_{x},{\varphi }_{x+1},{\varphi }_{x+2},{\varphi }_{x+3})$ and (24), (25) would result in some exactly solvable nonlinear functional equations. For instance, intriguing numerical results in the integrable lattice Landau–Lifshitz classical spin chain model suggest [10] existence of a nontrivial nonequilibrium phase transition from ballistic to diffusive steady-state and a nontrivial quasilocal conservation law in the ballistic regime, which, according to the results presented here, might be analytically treatable with appropriate integrable Markovian boundary baths.

For general classical integrable systems which are canonically defined via the zero-curvature condition in terms of a Lax pair, one should explore the possibility of connection between the patch tensors $L,X,R$ and the Lax operators generating equations of motion in the bulk. It should be noted, however, that compatibility/integrability condition between the deterministic bulk and stochastic boundaries would require all the patch tensors to explicitly depend on the Markov rates at the boundaries, just like in the model solved here. The question whether this can be accommodated for in terms of a single family of Lax operators, in a similar way as encoding the boundary dissipation in the complex auxiliary spin of the ${U}_{q}({{\mathfrak{sl}}}_{2})$ Lax operator in the case of boundary driven quantum XXZ chains [3, 9], remains an open problem for future research.

TP thanks E Ilievski for discussions and useful remarks on the manuscript. The work has been supported by the research grants P1-0044, J1-5439 and N1-0025 of Slovenian Research Agency (TP), and Spanish MICINN grants MTM2012-39101-C02-01 and MTM2015-63914-P (CM-M).