Rule 54: exactly solvable model of nonequilibrium statistical mechanics

Establishing the existence of a unique NESS and relaxation towards it from an arbitrary initial probability state vector amounts to showing the following statement [21]:

Theorem 1. The 2ⁿ × 2ⁿ matrix $\mathbb{U}$, equation (7), is irreducible and aperiodic for generic values of driving parameters, more precisely, for an open set 0 < α, β, γ, δ < 1 for conditional driving (15), and 0 < ζ, α, β, η, γ, δ < 1 for Bernoulli driving (13).

Proof. According to the Perron–Frobenius theorem [40], a finite, non-negative matrix $\mathbb{U}$ is irreducible, if for any pair of configurations $\underline{s},{\underline{s}}^{\prime }\in {\mathcal{C}}_{n}$ there exists a natural number ${t}_{0}\in \mathbb{N}$ such that the corresponding matrix element of the t₀th power is nonzero,

$$\begin{equation}{\left[{\mathbb{U}}^{{t}_{0}}\right]}_{{\underline{s}}^{\prime },\underline{s}}{ >}0.\end{equation} \tag{ 16 }$$

An irreducible matrix $\mathbb{U}$ is aperiodic if for any configuration $\underline{s}\in {\mathcal{C}}_{n}$, the greatest common divisor of the set of recurrence times {t_j } is 1,

$$\begin{equation}\mathrm{gcd}\left(\left\{{t}_{j}\in \mathbb{N};\enspace {\left[{\mathbb{U}}^{{t}_{j}}\right]}_{\underline{s},\underline{s}}{ >}0\right\}\right)=1.\end{equation} \tag{ 17 }$$

Let us first show irreducibility. By the definition of boundary propagators (cf (13) and (15)), the Markov matrix $\mathbb{U}$ connects each configuration $\underline{s}$ to exactly 4 other configurations ${\underline{s}}^{\prime }$, which exhibit all the possible configurations of boundary bits, $({s}_{1}^{\prime },{s}_{n}^{\prime })\in \left\{(0,0),(0,1),(1,0),(1,1)\right\}$. This holds for all values of parameters α, β, γ, δ (and ζ, η), except for the marginal case when some of the parameters are equal to 0 or 1, but this option is excluded by the assumption of the theorem. Let us now take a sufficiently large positive integer t₀, 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{equation}\underline{s}(0)\to \underline{s}(1)\to \underline{s}(2)\to \cdots \to \underline{s}({t}_{0}),\end{equation} \tag{ 18 }$$

which can be understood as a path through the Markov graph defined by positive elements of $\mathbb{U}$ that connects $\underline{s}$ and ${\underline{s}}^{\prime }$ in t₀ steps and implies ${\left[{\mathbb{U}}^{{t}_{0}}\right]}_{{\underline{s}}^{\prime },\underline{s}}{ >}0$ (see figure 4 for a 'self-contained' graphic illustration of the idea of proof). We are still free to choose the values of the boundary cells s_1,n (t) along the walk t ∈ {1, 2, ..., t₀ − 1} apart from the ends. For the first part of the walk t = 1, 2...t₊, up to some t₊ < t₀, we are fixing them with the rule

$$\begin{equation}{s}_{1}(t)={s}_{2}(t-1),\hspace{25.0pt}{s}_{n}(t)={s}_{n-1}(t-1).\end{equation} \tag{ 19 }$$

The evolution of the interior values of the cells s_x (t) for 1 < x < n and t ⩽ t₊ is then completely specified by the deterministic RCA54, while (19) provide the causal absorbing boundary conditions. Indeed, each time the boundary cell, say x = 1, gets occupied, s₁(t) = 1, the soliton is absorbed (see figure 4). As the solitons only move ballistically (at speed 1) and scatter pairwise (with time-lag 1), it is clear that a finite time scale ${t}_{+}\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\cdot \cdot \cdot ,0)$.

For the rest of the walk t ∈ {t₊ + 1, ...t₀} we need to show that alternative boundary rules exists, which create the configuration ${\underline{s}}^{\prime }$ out of the vacuum in another t₋ = t₀ − 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 (19) when the time runs backwards,

$$\begin{equation}{s}_{1}(t)={s}_{2}(t+1),\hspace{25.0pt}{s}_{n}(t)={s}_{n-1}(t+1),\quad \mathrm{f}\mathrm{o}\mathrm{r}\;t={t}_{0}-1,\cdot \cdot \cdot ,{t}_{0}-{t}_{-}.\end{equation} \tag{ 20 }$$

The entire walk then connects $\underline{s}$ to ${\underline{s}}^{\prime }$ in t₀ = t₊ + t₋ steps and implies ${\left[{\mathbb{U}}^{{t}_{0}}\right]}_{{\underline{s}}^{\prime },\underline{s}}{ >}0$, for an arbitrary pair $\underline{s},{\underline{s}}^{\prime }\in {\mathcal{C}}_{n}$, with the minimal possible integer t₀ in general depending 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 ${\left[{\mathbb{U}}^{{t}_{0}}\right]}_{\underline{s},\underline{s}}{ >}0$ for some t₀ = t₊ + t₋ depending on $\underline{s}$. However, in between annihilating the configuration $\underline{s}$ in t₊ time steps and then creating it again in another t₋ steps ⁴ , we can stay in the vacuum state for an arbitrary additional number of steps t_gap ⩾ 0. In this way the walk is increased by a segment of t_gap intermediate vacuum configurations, and we still have a non-zero matrix element,

$$\begin{equation}{\left[{\mathbb{U}}^{{t}_{0}+{t}_{\text{gap}}}\right]}_{\underline{s},\underline{s}}{ >}0.\end{equation} \tag{ 21 }$$

The greatest common divisor of the set $\left\{{t}_{0}+{t}_{\text{gap}} ;\enspace {t}_{\text{gap}}\in {\mathbb{Z}}_{+}\right\}$ is clearly 1, so we have shown aperiodicity. □

Zoom In Zoom Out Reset image size

In conclusion, the Perron–Frobenius theorem [40] guarantees that the NESS probability state vector p₀, satisfying the fixed point condition $\mathbb{U}{\mathbf{p}}_{0}={\mathbf{p}}_{0}$, is unique and all other eigenvalues of $\mathbb{U}$ lie strictly inside the unit circle. As a consequence, the Markov dynamics (6) is ergodic and mixing and an arbitrary initial probability state vector p(0) converges to NESS exponentially fast in t.

The NESS probability state vector p = p₀ can be split into even and odd time slice, satisfying a pair of fixed point equations

$$\begin{equation}{\mathbf{p}}^{\prime }={\mathbb{U}}_{\mathrm{e}}\mathbf{p},\hspace{25.0pt}\mathbf{p}={\mathbb{U}}_{\mathrm{o}}{\mathbf{p}}^{\prime }.\end{equation} \tag{ 22 }$$

We shall now explicitly construct the probability state vectors p and p', by solving equation (22) in terms of a simple ansatz, which we term a patch state ansatz (PSA) (illustrated in figure 5).

Zoom In Zoom Out Reset image size

For either of the two boundary driving families (13) and (15), the NESS solution $\mathbf{p},{\mathbf{p}}^{\prime }\in {\mathbb{R}}^{{2}^{n}}$ of the fixed point condition (22) can be written in the form

$$\begin{equation}\begin{aligned}{p}_{{s}_{1},{s}_{2},\cdot \cdot \cdot ,{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}}\cdot \cdot \cdot {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},\cdot \cdot \cdot ,{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 }\cdot \cdot \cdot {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{aligned}\end{equation} \tag{ 23 }$$

for some rank-4 and rank-3 tensors of strictly positive components X_ss'uu', ${X}_{s{s}^{\prime }u{u}^{\prime }}^{\prime }$, L_suu', ${L}_{su{u}^{\prime }}^{\prime }$, R_ss'u , ${R}_{s{s}^{\prime }u}^{\prime }$, with binary indices s, s', u, u' ∈ {0, 1}. To uniquely determine the algebraic expressions for these tensors in terms of the parameters of the model, one has to plug the ansatz (23) into the NESS condition (22).

First we find a minimal sufficient set of equations that the tensors have to satisfy, by fixing the normalization and taking into account the gauge symmetry of the ansatz. The normalization of the PSA (23) can be chosen such that

$$\begin{equation}{X}_{0000}={X}_{0000}^{\prime }=1,\hspace{25.0pt}{L}_{000}={R}_{000}^{\prime }=1.\end{equation} \tag{ 24 }$$

Clearly ${X}_{0000}={X}_{0000}^{\prime }$, otherwise the probabilities of the vacuum configurations p_00...0 and ${p}_{00\cdot \cdot \cdot 0}^{\prime }$ would scale differently with n which is not possible since they are directly connected by both even and odd propagators, i.e. ${\left[{\mathbb{U}}_{\mathrm{o}}\right]}_{(00\cdot \cdot \cdot 0)(00\cdot \cdot \cdot 0)},{\left[{\mathbb{U}}_{\mathrm{e}}\right]}_{(00\cdot \cdot \cdot 0)(00\cdot \cdot \cdot 0)}{ >}0$.

Let us now assume the ansatz (23) and write all components of equation (22), ${\left[{\mathbb{U}}_{\mathrm{e}}\mathbf{p}-{\mathbf{p}}^{\prime }\right]}_{\underline{s}}={\left[{\mathbb{U}}_{\mathrm{o}}{\mathbf{p}}^{\prime }-\mathbf{p}\right]}_{\underline{s}}=0$, pertaining to four-cluster configurations in the bulk of the form ⁵ $\underline{s}=({0}^{\left\{2k+1\right\}},s,{s}^{\prime },u,{u}^{\prime },{0}^{\left\{n-5-2k\right\}})$, for k = 0, 1, ..., m − 3, which results in the following two bulk equations:

$$\begin{equation}\begin{aligned}& {X}_{00s{s}^{\prime }}^{\prime }{X}_{s{s}^{\prime }u{u}^{\prime }}^{\prime }{X}_{u{u}^{\prime }00}^{\prime }{X}_{0000}^{\prime }={X}_{00\chi (0s{s}^{\prime }){s}^{\prime }}{X}_{\chi (0s{s}^{\prime }){s}^{\prime }\chi ({s}^{\prime }u{u}^{\prime }){u}^{\prime }}{X}_{\chi ({s}^{\prime }u{u}^{\prime }){u}^{\prime }\chi ({u}^{\prime }00)0}{X}_{\chi ({u}^{\prime }00)000},\\ & {X}_{0000}{X}_{00s{s}^{\prime }}{X}_{s{s}^{\prime }u{u}^{\prime }}{X}_{u{u}^{\prime }00}={X}_{000\chi (00s)}^{\prime }{X}_{0\chi (00s)s\chi (s{s}^{\prime }u)}^{\prime }{X}_{s\chi (s{s}^{\prime }u)u\chi (u{u}^{\prime }0)}^{\prime }{X}_{u\chi (u{u}^{\prime }0)00}^{\prime }.\end{aligned}\end{equation} \tag{ 25 }$$

Similarly, considering three-cluster configurations near each boundary, $\underline{s}=({v}^{\prime },s,{s}^{\prime },{0}^{\left\{n-3\right\}})$ and $\underline{s}=({0}^{\left\{n-3\right\}},s,{s}^{\prime },u)$, we obtain a set of four boundary equations,

$$\begin{equation}\begin{aligned}{L}_{{v}^{\prime }s{s}^{\prime }}^{\prime }{X}_{s{s}^{\prime }00}^{\prime }{X}_{0000}^{\prime }{R}_{000}^{\prime }& ={L}_{{v}^{\prime }\chi ({s}^{\prime }s{s}^{\prime }){s}^{\prime }}{X}_{\chi ({v}^{\prime }s{s}^{\prime }){s}^{\prime }\chi ({s}^{\prime }00)0}{X}_{\chi ({s}^{\prime }00)000}\frac{{R}_{000}+{R}_{001}}{2},\\ {L}_{000}^{\prime }{X}_{00s{s}^{\prime }}^{\prime }{R}_{s{s}^{\prime }u}^{\prime }& ={L}_{000}\sum\limits _{{t}^{\prime },t}{P}_{({s}^{\prime },u),({t}^{\prime },t)}^{\mathrm{R}}{X}_{00\chi (0s{t}^{\prime }){t}^{\prime }}{R}_{\chi (0s{t}^{\prime }){t}^{\prime }t},\\ {L}_{{v}^{\prime }s{s}^{\prime }}{X}_{s{s}^{\prime }00}{R}_{000}& =\sum\limits _{{t}^{\prime },t}{P}_{({v}^{\prime },s),({t}^{\prime },t)}^{\mathrm{L}}{L}_{{t}^{\prime }t\chi (t{s}^{\prime }0)}^{\prime }{X}_{t\chi (t{s}^{\prime }0)00}^{\prime }{R}_{000}^{\prime },\\ {L}_{000}{X}_{0000}{X}_{00s{s}^{\prime }}{R}_{s{s}^{\prime }u}& =\frac{{L}_{000}^{\prime }+{L}_{100}^{\prime }}{2}{X}_{000\chi (00s)}^{\prime }{X}_{0\chi (00s)s\chi (s{s}^{\prime }u)}^{\prime }{R}_{s\chi (s{s}^{\prime }u)u}^{\prime }.\end{aligned}\end{equation} \tag{ 26 }$$

The total number of 2 × 16 + 4 × 8 − 4 = 60 unknowns can be further reduced by exploiting the following gauge symmetry

$$\begin{equation}\begin{cases}{X}_{s{s}^{\prime }t{t}^{\prime }}\to {f}_{s{s}^{\prime }}{X}_{s{s}^{\prime }t{t}^{\prime }}{f}_{t{t}^{\prime }}^{-1},\\ \enspace {L}_{{s}^{\prime }t{t}^{\prime }}\to {L}_{{s}^{\prime }t{t}^{\prime }}{f}_{t{t}^{\prime }}^{-1},\\ \enspace {R}_{s{s}^{\prime }t}\to {f}_{s{s}^{\prime }}{R}_{s{s}^{\prime }t},\end{cases}\hspace{15.0pt}\begin{cases}{X}_{s{s}^{\prime }t{t}^{\prime }}^{\prime }\to {g}_{s{s}^{\prime }}{X}_{s{s}^{\prime }t{t}^{\prime }}^{\prime }{g}_{t{t}^{\prime }}^{-1},\\ \enspace {L}_{{s}^{\prime }t{t}^{\prime }}^{\prime }\to {L}_{{s}^{\prime }t{t}^{\prime }}^{\prime }{g}_{t{t}^{\prime }}^{-1},\\ \enspace {R}_{s{s}^{\prime }t}^{\prime }\to {g}_{s{s}^{\prime }}{R}_{s{s}^{\prime }t}^{\prime },\end{cases}\end{equation} \tag{ 27 }$$

which conserves the patch ansatz (23), as well as the defining equations (25) and (26) for arbitrary nonzero gauge 'fields' f_ss', g_ss'.

While a detailed analysis and explicit expressions for the PSA tensors in terms of boundary parameters can be found in [21, 22], here we only discuss a generic solution form of the bulk part of the equations. Specifically, the bulk equations (25) can be solved independently (before incorporating boundary conditions), and the solution can be parametrised uniquely—up to a choice of gauge (27)—in terms of two free (spectral) parameters ξ, ω

$$\begin{equation}{X}_{s{s}^{\prime }t{t}^{\prime }}={T}_{(s{s}^{\prime }),(t{t}^{\prime })}(\xi ,\omega ),\hspace{10.0pt}{X}_{s{s}^{\prime }t{t}^{\prime }}^{\prime }={T}_{(s{s}^{\prime }),(t{t}^{\prime })}(\omega ,\xi ),\hspace{10.0pt}T(\xi ,\omega ) = \left[\begin{matrix}\hfill 1\hfill & \hfill 1\hfill & \hfill \xi \hfill & \hfill 1\hfill \\ \hfill \xi \omega \hfill & \hfill \xi \omega \hfill & \hfill 1\hfill & \hfill \omega \hfill \\ \hfill \omega \hfill & \hfill \omega \hfill & \hfill \xi \omega \hfill & \hfill \omega \hfill \\ \hfill \xi \hfill & \hfill \xi \hfill & \hfill \xi \hfill & \hfill \xi \omega \hfill \end{matrix}\right] .\end{equation} \tag{ 28 }$$

The boundary equations (26) then fix the spectral parameters ξ, ω as functions of the boundary driving parameters α, β, γ, δ (and ξ, η) (see section 4). The remarkable fact is that the solutions do not explicitly depend on the system size n, hence it is clear that such an NESS can only describe ballistic transport (i.e. net soliton current independent of n).

The 4 × 4 matrix T(ξ, ω), defined in equation (28), can be interpreted as a two-parametric transfer matrix generating the local charges of the RCA54 model. To simplify the discussion, let us assume periodic boundary conditions, n + 1 ≡ 1, and consider a steady state vector p(ξ, ω) whose components are given in terms of the transfer matrix as,

$$\begin{equation}{p}_{{s}_{1},{s}_{2},\cdot \cdot \cdot ,{s}_{n}}(\xi ,\omega )={T}_{({s}_{1},{s}_{2}),({s}_{3},{s}_{4})}(\xi ,\omega ){T}_{({s}_{3},{s}_{4}),({s}_{5},{s}_{6})}(\xi ,\omega )\cdot \cdot \cdot {T}_{({s}_{n-1},{s}_{n}),({s}_{1},{s}_{2})}(\xi ,\omega ).\end{equation} \tag{ 29 }$$

For any pair of positive real parameters ξ, ω > 0, p(ξ, ω) represents the statistical ensemble—state that is invariant under time translation (as well as translationally invariant in space)

$$\begin{equation}{\mathbb{U}}_{\mathrm{e}}\mathbf{p}(\xi ,\omega )=\mathbf{p}(\omega ,\xi ),\quad {\mathbb{U}}_{\mathrm{o}}\mathbf{p}(\omega ,\xi )=\mathbf{p}(\xi ,\omega ),\quad \mathbb{U}\mathbf{p}(\xi ,\omega )=\mathbf{p}(\xi ,\omega ),\end{equation} \tag{ 30 }$$

whose partition sum is given simply in terms of powers of the transfer matrix as

$$\begin{equation}\sum\limits _{\underline{s}}{p}_{\underline{s}}(\xi ,\omega )=\enspace \text{tr}\enspace {[T(\xi ,\omega )]}^{n/2}.\end{equation} \tag{ 31 }$$

Parameters log ξ, log ω can be interpreted as chemical potentials corresponding to the conserved numbers of left and right movers, N₋ and N₊,

$$\begin{equation}{p}_{\underline{s}}(\xi ,\omega )\propto {\xi }^{{N}_{-}(\underline{s})}{\omega }^{{N}_{+}(\underline{s})}.\end{equation} \tag{ 32 }$$

Specifically, the latter can be generated directly from the logarithmic derivatives of the transfer matrix

$$\begin{equation}\begin{aligned}{N}_{+}(\underline{s})=\frac{\mathrm{d}}{\mathrm{d}\enspace \mathrm{log}\enspace \xi }\enspace \mathrm{log}\enspace \prod\limits _{j=1}^{n/2}{\left.{T}_{({s}_{2j-1},{s}_{2j}),({s}_{2j+1},{s}_{2j+2})}(\xi ,\omega )\right\vert }_{\xi =\omega =1},\\ {N}_{-}(\underline{s})=\frac{\mathrm{d}}{\mathrm{d}\enspace \mathrm{log}\enspace \omega }\enspace \mathrm{log}\enspace \prod\limits _{j=1}^{n/2}{\left.{T}_{({s}_{2j-1},{s}_{2j}),({s}_{2j+1},{s}_{2j+2})}(\xi ,\omega )\right\vert }_{\xi =\omega =1},\end{aligned}\end{equation} \tag{ 33 }$$

which can clearly be written as extensive sums of local densities

$$\begin{equation}{N}_{{\pm}}=\sum\limits _{j=1}^{n/2}{n}_{{\pm}}^{j},\hspace{25.0pt}\begin{aligned}{n}_{+}^{j}(\underline{s})={\left.\frac{\mathrm{d}}{\mathrm{d}\xi }{T}_{({s}_{2j-1},{s}_{2j}),({s}_{2j+1},{s}_{2j+2})}\right\vert }_{\xi =\omega =1},\\ {n}_{-}^{j}(\underline{s})={\left.\frac{\mathrm{d}}{\mathrm{d}\omega }{T}_{({s}_{2j-1},{s}_{2j}),({s}_{2j+1},{s}_{2j+2})}\right\vert }_{\xi =\omega =1}.\end{aligned}\end{equation} \tag{ 34 }$$

These are the elementary local charges of RCA54 and will be discussed later in the context of hydrodynamics of the model. We note that N_± by no means represent a complete set of local charges. Specifically, searching for all local charges with densities of support size ℓ > 3, we find of the order of Fibonacci number F_ℓ = F_ℓ−1 + F_ℓ−2 of trivially conserved charges corresponding to F_ℓ non-interacting unidirectional soliton configurations. See the discussion in [41] for the details.

Following [23, 25] we take a staggered MPA for the NESS in the form,

$$\begin{equation}\mathbf{p}=\left\langle {\mathbf{l}}_{1}\right\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }{\mathbf{W}}_{4}{\mathbf{W}}_{5}^{\prime }\cdot \cdot \cdot {\mathbf{W}}_{n-3}^{\prime }{\mathbf{W}}_{n-2}\vert {\mathbf{r}}_{n-1,n}\rangle ,\end{equation} \tag{ 38 }$$

$$\begin{equation}{\mathbf{p}}^{\prime }=\left\langle {\mathbf{l}}_{12}^{\prime }\right\vert {\mathbf{W}}_{3}{\mathbf{W}}_{4}^{\prime }{\mathbf{W}}_{5}{\mathbf{W}}_{6}^{\prime }\cdot \cdot \cdot {\mathbf{W}}_{n-2}^{\prime }{\mathbf{W}}_{n-1}\vert {\mathbf{r}}_{n}^{\prime }\rangle .\end{equation} \tag{ 39 }$$

We will use this as an ansatz to solve the split eigenvalue equation (36). We require that the operators ${\mathit{W}}_{s},{\mathit{W}}_{s}^{\prime }$ and some operator S satisfy the following bulk algebraic relation,

$$\begin{equation}{U}_{123}{\mathbf{W}}_{1}S\enspace {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }={\mathbf{W}}_{1}{\mathbf{W}}_{2}^{\prime }{\mathbf{W}}_{3}S.\end{equation} \tag{ 40 }$$

Component-wise the above vector equation reads

$$\begin{equation}{W}_{s}S{W}_{\chi (s,{s}^{\prime },{s}^{{\prime\prime}})}{W}_{{s}^{{\prime\prime}}}^{\prime }={W}_{s}{W}_{{s}^{\prime }}^{\prime }{W}_{{s}^{{\prime\prime}}}S,\quad s,{s}^{\prime },{s}^{{\prime\prime}}\in \left\{0,1\right\}.\end{equation} \tag{ 41 }$$

A representation of this algebra can be found on a three-dimensional auxiliary vector space ${\mathcal{V}}_{a}={\mathbb{C}}^{3}$ in terms of matrices

$$\begin{equation}{W}_{0}(\xi ,\omega )=\left[\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \xi \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}{W}_{1}(\xi ,\omega )=\left[\begin{matrix}\hfill 0\hfill & \hfill \xi \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \omega \hfill \end{matrix}\right],\hspace{25.0pt}S=\left[\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right],\end{equation} \tag{ 42 }$$

while the other set of auxiliary space matrices ${W}_{s}^{\prime }(\xi ,\omega )$ is obtained from W_s (ξ, ω) by swapping the parameters,

$$\begin{equation}{W}_{s}^{\prime }(\xi ,\omega )={W}_{s}(\omega ,\xi ),\quad s\in \left\{0,1\right\}.\end{equation} \tag{ 43 }$$

The representation is parametrised in terms of two (so far) free variables ξ and ω, which (as we argue later) correspond precisely to the variables parametrising the transfer matrix introduced in (28), therefore we shall refer to them as spectral parameters. We note that ${S}^{2}=\mathbb{1}$, which together with ${U}^{2}=\mathbb{1}$ and the bijection between the two sets of matrices implies a dual bulk relation,

$$\begin{equation}{U}_{123}{\mathbf{W}}_{1}^{\prime }{\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }S={\mathbf{W}}_{1}^{\prime }S\enspace {\mathbf{W}}_{2}^{\prime }{\mathbf{W}}_{3}.\end{equation} \tag{ 44 }$$

We now turn to the boundaries. We demand that the spectral parameters ξ, ω, eigenvalue parameters Λ_L, Λ_R and the boundary vectors $\left\langle {\mathbf{l}}_{1}\right\vert$, $\left\langle {\mathbf{l}}_{12}^{\prime }\right\vert$, |r₁₂⟩, $\vert {\mathbf{r}}_{1}^{\prime }\rangle$, satisfy the following set of equations,

$$\begin{equation}{U}_{12}^{\mathrm{L}}\left\langle {\mathbf{l}}_{12}^{\prime }\right\vert ={{\Lambda}}_{\mathrm{L}}\left\langle {\mathbf{l}}_{1}\right\vert {\mathbf{W}}_{2}S,\end{equation} \tag{ 45 }$$

$$\begin{equation}{U}_{123}\left\langle {\mathbf{l}}_{1}\right\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }=\left\langle {\mathbf{l}}_{12}^{\prime }\right\vert {\mathbf{W}}_{3}S,\end{equation} \tag{ 46 }$$

$$\begin{equation}{U}_{12}^{\mathrm{R}}\vert {\mathbf{r}}_{12}\rangle ={{\Lambda}}_{\mathrm{R}}{\mathbf{W}}_{1}^{\prime }S\vert {\mathbf{r}}_{2}^{\prime }\rangle ,\end{equation} \tag{ 47 }$$

$$\begin{equation}{U}_{123}{\mathbf{W}}_{1}^{\prime }{\mathbf{W}}_{2}\vert {\mathbf{r}}_{3}^{\prime }\rangle ={\mathbf{W}}_{1}^{\prime }S\vert {\mathbf{r}}_{23}\rangle .\end{equation} \tag{ 48 }$$

If these equations are satisfied, the pair of states p, p' written in terms of MPAs (38) and (39), solves the staggered eigenvalue equations (36). This can be straightforwardly verified by plugging the MPAs into e.g. the first equation of the set (36), and apply the appropriate relations to transform p into p',

$$\begin{align}{\mathbb{U}}_{\mathrm{e}}\mathbf{p}& ={U}_{123}\cdot \cdot \cdot {U}_{n-5\enspace n-4\enspace n-3}{U}_{n-3\enspace n-2\enspace n-1}{U}_{n-1\enspace n}^{\mathrm{R}}\left\langle {\mathbf{l}}_{1}\right\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }\cdot \cdot \cdot {\mathbf{W}}_{n-4}{\mathbf{W}}_{n-3}^{\prime }{\mathbf{W}}_{n-2}\left\vert {\mathbf{r}}_{n-1\enspace n}\right\rangle \\ & ={{\Lambda}}_{\mathrm{R}}{U}_{123}\cdot \cdot \cdot {U}_{n-5\enspace n-4\enspace n-3}{U}_{n-3\enspace n-2\enspace n-1}\left\langle {\mathbf{l}}_{1}\left\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }\cdot \cdot \cdot {\mathbf{W}}_{n-4}{\mathbf{W}}_{n-3}^{\prime }{\mathbf{W}}_{n-2}{\mathbf{W}}_{n-1}^{\prime }S\right\vert {\mathbf{r}}_{n}^{\prime }\right\rangle \\ & ={{\Lambda}}_{\mathrm{R}}{U}_{123}\cdot \cdot \cdot {U}_{n-5\enspace n-4\enspace n-3}\left\langle {\mathbf{l}}_{1}\left\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }\cdot \cdot \cdot {\mathbf{W}}_{n-4}{\mathbf{W}}_{n-3}^{\prime }S{\mathbf{W}}_{n-2}^{\prime }{\mathbf{W}}_{n-1}\right\vert {\mathbf{r}}_{n}^{\prime }\right\rangle \\ & =\cdot \cdot \cdot ={{\Lambda}}_{\mathrm{R}}{U}_{123}\left\langle {\mathbf{l}}_{1}\left\vert {\mathbf{W}}_{2}{\mathbf{W}}_{3}^{\prime }S{\mathbf{W}}_{4}^{\prime }{\mathbf{W}}_{5}\cdot \cdot \cdot {\mathbf{W}}_{n-1}\right\vert {\mathbf{r}}_{n}^{\prime }\right\rangle \\ & ={{\Lambda}}_{\mathrm{R}}\left\langle {\mathbf{l}}_{12}^{\prime }\left\vert {\mathbf{W}}_{3}{\mathbf{W}}_{4}^{\prime }{\mathbf{W}}_{5}\cdot \cdot \cdot {\mathbf{W}}_{n-1}\right\vert {\mathbf{r}}_{n}^{\prime }\right\rangle ={{\Lambda}}_{\mathrm{R}}{\mathbf{p}}^{\prime }.\end{align} \tag{ 49 }$$

To get to the second line, we apply the boundary propagator ${U}_{n-1\enspace n}^{\mathrm{R}}$, which, according to relation (47), replaces the right boundary vector $\left\vert {\mathbf{r}}_{n-1\enspace n}^{\prime }\right\rangle$ with $\left\vert {\mathbf{r}}_{n}\right\rangle$, and at the same time introduces the delimiter matrix S. We continue by applying the dual bulk relation (44), which moves S to the left, while the roles of ${\mathbf{W}}_{s}^{\prime }$ and W_s are exchanged (fourth line). We finish by absorbing the matrix S at the left and finally change $\left\langle {\mathbf{l}}_{1}\right\vert$ into $\left\langle {\mathbf{l}}_{12}^{\prime }\right\vert$. The second part of (36) is proven analogously.

Solving the left boundary equations (45) and (46) for ξ, ω and boundary vectors $\langle {l}_{s}\vert ,\langle {l}_{s{s}^{\prime }}^{\prime }\vert$ (given explicitly in appendix A) we obtain

$$\begin{equation}\xi =\frac{(\alpha +\beta -1)-\beta {{\Lambda}}_{\mathrm{L}}}{(\beta -1){{\Lambda}}_{\mathrm{L}}^{2}},\hspace{25.0pt}\omega =\frac{{{\Lambda}}_{\mathrm{L}}(\alpha -{{\Lambda}}_{\mathrm{L}})}{(\beta -1)}.\end{equation} \tag{ 50 }$$

Similarly, solving the right boundary equations (47) and (48) for ξ, ω and the $\vert {r}_{s}^{\prime }\rangle ,\vert {r}_{s{s}^{\prime }}\rangle$ (again listed in appendix A) we get

$$\begin{equation}\xi =\frac{{{\Lambda}}_{\mathrm{R}}(\gamma -{{\Lambda}}_{\mathrm{R}})}{(\delta -1)},\hspace{25.0pt}\omega =\frac{(\gamma +\delta -1)-\delta {{\Lambda}}_{\mathrm{R}}}{(\delta -1){{\Lambda}}_{\mathrm{R}}^{2}}.\end{equation} \tag{ 51 }$$

These and all explicit result below are written for the specific case of conditional driving (15), while results for the Bernoulli family are obtained fully analogously [23]. Equating ξ and ω from these pairs of equations yields an eigenvalue equation, which expressed in terms of the eigenvalue Λ = Λ_LΛ_R reads

$$\begin{equation}({\Lambda}-1)\left({{\Lambda}}^{3}+(1-\alpha \gamma ){{\Lambda}}^{2}+\beta \delta \enspace {\Lambda}-(\alpha +\beta -1)(\gamma +\delta -1)\left({\Lambda}+1\right)\right)=0.\end{equation} \tag{ 52 }$$

Importantly, Λ = 1 is always a solution. The corresponding eigenvector is the NESS. The remainder is a cubic polynomial giving three other eigenvalues corresponding to three decay modes. These eigenvalues do not depend on the system size n. We refer to these four eigenvalues as the NESS orbital.

In order to generalise the results from the previous subsection to other eigenvalues and eigenvectors (decay modes) of the Markov matrix $\mathbb{U}$ we will follow the spirit of the standard coordinate Bethe ansatz (see e.g. [46]) for finding solutions to eigenstates containing interacting particles. However, we need to adapt it to a coordinate MPA picture similar to the one used for the asymmetric simple exclusion processes [47]. As we will see the leading decay modes can be understood in terms of localized particle excitations of the NESS, the latter being analogous to the vacuum state. The leading decay modes will thus be a superposition of single-particle excitations.

The leading decay mode is the eigenvector of $\mathbb{U}$ corresponding to the eigenvalue Λ₁ with the real part closest to 1,

$$\begin{equation}\underset{j{\geqslant}1}{\mathrm{min}}\left(1-\left\vert \mathfrak{R}({{\Lambda}}_{j})\right\vert \right)=1-\left\vert \mathfrak{R}({{\Lambda}}_{1})\right\vert .\end{equation} \tag{ 53 }$$

Apart from the eigenvectors with eigenvalues on the unit circle, it is the most dynamically relevant eigenvector, since it dominates the long-time asymptotic relaxation. In this subsection, we will provide a compact MPA representation of the eigenvectors corresponding to the orbital that contains the leading decay mode. The generalisation to other orbitals remains an open question, although the complete set of eigenvalues can be predicted by a simple conjecture.

The starting point is to double the auxiliary space, which now consists of two copies of the original auxiliary space ${\mathcal{V}}_{a}^{\prime }={\mathcal{V}}_{a}\oplus {\mathcal{V}}_{a}={\mathbb{C}}^{2}\otimes {\mathcal{V}}_{a}$. We generalise the W operator to the doubled auxiliary space ${\mathcal{V}}_{a}^{\prime }$ as,

$$\begin{equation}\begin{aligned}{\hat{W}}_{s}={e}_{11}\otimes {W}_{s}(\xi z,\omega /z)+{e}_{22}\otimes {W}_{s}(\xi /z,\omega z),\\ {\hat{W}}_{s}^{\prime }={e}_{11}\otimes {W}_{s}^{\prime }(\xi z,\omega /z)+{e}_{22}\otimes {W}_{s}^{\prime }(\xi /z,\omega z),\end{aligned}\end{equation} \tag{ 54 }$$

where ${e}_{ij}=\left\vert i\right\rangle \left\langle j\right\vert$, i, j ∈ {1, 2} is the Weyl basis of 2 × 2 matrices. We define a set of operators ${\hat{F}}^{(k)}$, ${\hat{F}}^{\prime (k)}$, ${\hat{G}}^{(k)}$, ${\hat{G}}^{\prime (k)}$, ${\hat{K}}^{(k)}$, ${\hat{L}}^{\prime (k)}$, acting on ${\mathcal{V}}_{a}^{\prime }$,

$$\begin{align}{\hat{F}}^{(k)}& =\mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{F}_{+}+{c}_{-}{z}^{-k}{F}_{-}}{\xi \omega -1},\hspace{25.0pt}{\hat{F}}^{\prime (k)}=\mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{F}_{+}^{\prime }+{c}_{-}{z}^{-k}{F}_{-}^{\prime }}{\xi \omega -1},\\ \hspace{25.0pt}{\hat{G}}^{(k)}& =\mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{G}_{+}+{c}_{-}{z}^{-k}{G}_{-}}{\xi \omega -1},\hspace{22.0pt}{\hat{G}}^{\prime (k)}=\mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{G}_{+}^{\prime }+{c}_{-}{z}^{-k}{G}_{-}^{\prime }}{\xi \omega -1},\\ {\hat{K}}^{(k)}& =\mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{K}_{+}+{c}_{-}{z}^{-k}{K}_{-}}{\xi \omega -1},\\ {\hat{L}}^{(k)}& =(z{e}_{11}+{z}^{-1}{e}_{22})\otimes \mathbb{1}+{e}_{12}\otimes \frac{{c}_{+}{z}^{k}{L}_{+}+{c}_{-}{z}^{-k}{L}_{-}}{\xi \omega -1},\end{align} \tag{ 55 }$$

where we introduced the operators ${F}_{{\pm}}^{(\prime )}$, ${G}_{{\pm}}^{(\prime )}$, K_±, L_± that act on ${\mathcal{V}}_{a}$, i.e. only on a single copy of the original auxiliary space. The explicit representation of these operators is given in reference [23] in terms of a four dimensional auxiliary space ${\mathcal{V}}_{a}={\mathbb{C}}^{4}$. Note that unlike the MPA solution corresponding to the NESS orbital, we have so far been unable to find a three-dimensional representation of auxiliary space operators, so a four-dimensional single-copy auxiliary space is assumed here. The operators (55) linearly depend on two complex coefficients c₊, c₋, and in the inhomogeneous MPA (introduced below), the complex variable z plays the role of the multiplicative quasi-momentum of the quasi-particle excitation on top of the NESS, which are created by the operators (55).

We now introduce the following inhomogeneous (site-dependent) MPA

$$\begin{equation}\mathbf{p}=\left\langle {\hat{\mathbf{l}}}_{1}\right\vert {\hat{L}}^{(0)}{\hat{\mathbf{W}}}_{2}^{\prime }\hat{S}{\hat{F}}^{(1)}{\hat{\mathbf{W}}}_{3}^{\prime }{\hat{G}}^{(2)}{\hat{\mathbf{W}}}_{4}\cdot \cdot \cdot {\hat{F}}^{(n-5)}{\hat{\mathbf{W}}}_{n-3}^{\prime }{\hat{G}}^{(n-4)}{\hat{\mathbf{W}}}_{n-2}\left\vert {\hat{\mathbf{r}}}_{n-1,n}\right\rangle ,\end{equation} \tag{ 56 }$$

$$\begin{equation}{\mathbf{p}}^{\prime }=\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{F}}^{\prime (1)}{\hat{\mathbf{W}}}_{3}{\hat{G}}^{\prime (2)}{\hat{\mathbf{W}}}_{4}^{\prime }{\hat{F}}^{\prime (3)}{\hat{\mathbf{W}}}_{5}\cdot \cdot \cdot {\hat{G}}^{\prime (n-4)}{\hat{\mathbf{W}}}_{n-2}^{\prime }{\hat{K}}^{(n-3)}{\hat{\mathbf{W}}}_{n-1}\hat{S}\left\vert {\hat{\mathbf{r}}}_{n}^{\prime }\right\rangle ,\end{equation} \tag{ 57 }$$

and p is the eigenvector of $\mathbb{U}={\mathbb{U}}_{\mathrm{e}}{\mathbb{U}}_{\mathrm{o}}$ (8) and (9) with eigenvalue Λ = Λ_LΛ_R, like before. The interpretation of z as the momentum of a single excitation should now be clear. Due to the presence of the off-diagonal nilpotent e₁₂ in the site dependent matrices ${\hat{F}}^{(k)},{\hat{G}}^{(k)}$ these can only create at most one excitation on the NESS in the MPA (56) and (57). Furthermore, since the site dependent matrices contain both terms with z and z⁻¹, this MPA is like a standing wave. We deform the bulk relation and its dual, (40) and (44), to an inhomogeneous form—where a free integer k denotes a spatial coordinate—which our newly introduced operators have to satisfy

$$\begin{equation}\begin{aligned}{U}_{123}{\hat{K}}^{(k-1)}{\hat{\mathbf{W}}}_{1}\hat{S}{\hat{G}}^{(k)}{\hat{\mathbf{W}}}_{2}{\hat{F}}^{(k+1)}{\hat{\mathbf{W}}}_{3}^{\prime }={\hat{F}}^{\prime (k-1)}{\hat{\mathbf{W}}}_{1}{\hat{G}}^{\prime (k)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{K}}^{(k+1)}{\hat{\mathbf{W}}}_{3}\hat{S},\\ {U}_{123}{\hat{G}}^{\prime (k-1)}{\hat{\mathbf{W}}}_{1}^{\prime }{\hat{F}}^{\prime (k)}{\hat{\mathbf{W}}}_{2}{\hat{L}}^{(k+1)}{\hat{\mathbf{W}}}_{3}^{\prime }\hat{S}={\hat{L}}^{(k-1)}{\hat{\mathbf{W}}}_{1}^{\prime }\hat{S}{\hat{F}}^{(k)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{G}}^{(k+1)}{\hat{\mathbf{W}}}_{3},\end{aligned}\end{equation} \tag{ 58 }$$

where $\hat{S}=\mathbb{1}\otimes S$. Analogously as before we demand the boundary vectors and the other parameters to provide a nontrivial solution to a set of generalised boundary equations,

$$\begin{equation}{U}_{12}^{\mathrm{L}}\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert ={{\Lambda}}_{\mathrm{L}}\left\langle {\hat{\mathbf{l}}}_{1}\right\vert {\hat{G}}^{\prime (0)}{\hat{\mathbf{W}}}_{2}^{\prime },\end{equation} \tag{ 59 }$$

$$\begin{equation}{U}_{123}\left\langle {\hat{\mathbf{l}}}_{1}\right\vert {\hat{L}}^{(0)}{\hat{\mathbf{W}}}_{2}^{\prime }\hat{S}{\hat{F}}^{(1)}{\hat{\mathbf{W}}}_{3}^{\prime }=\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{K}}^{(1)}{\hat{\mathbf{W}}}_{3}\hat{S},\end{equation} \tag{ 60 }$$

$$\begin{equation}{U}_{12}^{\mathrm{R}}\left\vert {\hat{\mathbf{r}}}_{12}\right\rangle ={\hat{F}}^{(n-3)}{\hat{\mathbf{W}}}_{1}^{\prime }\hat{S}\left\vert {\hat{\mathbf{r}}}_{2}^{\prime }\right\rangle ,\end{equation} \tag{ 61 }$$

$$\begin{equation}{U}_{123}{\hat{G}}^{(n-4)}{\hat{\mathbf{W}}}_{1}^{\prime }{\hat{K}}^{(n-3)}{\hat{\mathbf{W}}}_{2}\hat{S}\left\vert {\hat{\mathbf{r}}}_{3}^{\prime }\right\rangle ={{\Lambda}}_{\mathrm{R}}{\hat{L}}^{(n-4)}{\hat{\mathbf{W}}}_{1}^{\prime }\hat{S}\left\vert {\hat{\mathbf{r}}}_{23}\right\rangle .\end{equation} \tag{ 62 }$$

To understand why p is indeed an eigenvector, let us for clarity consider an example with fixed n = 8 spin sites,

$$\begin{equation*}{\mathbb{U}}_{\mathrm{e}}\mathbf{p}={U}_{123}{U}_{345}{U}_{567}{U}_{78}^{\mathrm{R}}\left\langle {\hat{\mathbf{l}}}_{1}\right\vert {\hat{L}}^{(0)}{\hat{\mathbf{W}}}_{2}^{\prime }\hat{S}{\hat{F}}^{(1)}{\hat{\mathbf{W}}}_{3}^{\prime }{\hat{G}}^{(2)}{\hat{\mathbf{W}}}_{4}{\hat{F}}^{(3)}{\hat{\mathbf{W}}}_{5}^{\prime }{\hat{G}}^{(4)}{\hat{\mathbf{W}}}_{6}\left\vert {\hat{\mathbf{r}}}_{78}\right\rangle =\cdot \cdot \cdot \end{equation*}$$

We start by first acting with U₁₂₃ (note that all the local operators commute, cf (12), and can be applied in any order), which using (60) creates ${\hat{K}}^{(1)}{\hat{\mathbf{W}}}_{3}\hat{S}$ giving,

$$\begin{equation*}{\mathbb{U}}_{\mathrm{e}}\mathbf{p}={U}_{345}{U}_{567}{U}_{78}^{\mathrm{R}}\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{K}}^{(1)}{\hat{\mathbf{W}}}_{3}\hat{S}{\hat{G}}^{(2)}{\hat{\mathbf{W}}}_{4}{\hat{F}}^{(3)}{\hat{\mathbf{W}}}_{5}^{\prime }{\hat{G}}^{(4)}{\hat{\mathbf{W}}}_{6}\left\vert {\hat{\mathbf{r}}}_{78}\right\rangle .\end{equation*}$$

The subsequent U's in ${\mathbb{U}}_{\mathrm{e}}$ transfer ${\hat{K}}^{(1)}{\hat{\mathbf{W}}}_{3}\hat{S}$ to the penultimate right-hand side via the bulk relation (58),

$$\begin{equation*}\cdot \cdot \cdot ={U}_{567}{U}_{78}^{\mathrm{R}}\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{F}}^{\prime (1)}{\hat{\mathbf{W}}}_{1}{\hat{G}}^{\prime (2)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{K}}^{(3)}{\hat{\mathbf{W}}}_{3}\hat{S}{\hat{G}}^{(4)}{\hat{\mathbf{W}}}_{6}\left\vert {\hat{\mathbf{r}}}_{78}\right\rangle =\cdot \cdot \cdot \end{equation*}$$

Prior to acting with the final U_{n−3,n−2,n−1}, we apply U^R so it creates ${\hat{F}}^{(n-3)}{\hat{\mathbf{W}}}_{n-1}^{\prime }\hat{S}$ at the right end,

$$\begin{equation*}\cdot \cdot \cdot ={U}_{567}\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{F}}^{\prime (1)}{\hat{\mathbf{W}}}_{1}{\hat{G}}^{\prime (2)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{K}}^{(3)}{\hat{\mathbf{W}}}_{3}\hat{S}{\hat{G}}^{(4)}{\hat{\mathbf{W}}}_{6}{\hat{F}}^{(5)}{\hat{\mathbf{W}}}_{7}^{\prime }\hat{S}\left\vert {\hat{\mathbf{r}}}_{8}^{\prime }\right\rangle =\cdot \cdot \cdot \end{equation*}$$

We then use the top equation (58) for the final U_{n−3,n−2,n−1} to move ${\hat{K}}^{(n-5)}{\hat{\mathbf{W}}}_{n-3}\hat{S}$ to the far end into ${\hat{K}}^{(n-3)}{\hat{\mathbf{W}}}_{n-1}\hat{S}$ and we thus finally obtain the form (57) of p',

$$\begin{equation}\cdot \cdot \cdot =\left\langle {\hat{\mathbf{l}}}_{12}^{\prime }\right\vert {\hat{F}}^{\prime (1)}{\hat{\mathbf{W}}}_{1}{\hat{G}}^{\prime (2)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{F}}^{\prime (3)}{\hat{\mathbf{W}}}_{1}{\hat{G}}^{\prime (4)}{\hat{\mathbf{W}}}_{2}^{\prime }{\hat{K}}^{(5)}{\hat{\mathbf{W}}}_{3}{\hat{S}}^{2}\left\vert {\hat{\mathbf{r}}}_{8}^{\prime }\right\rangle ={\mathbf{p}}^{\prime }.\end{equation} \tag{ 63 }$$

The odd-part of the propagator ${\mathbb{U}}_{\mathrm{o}}$ acts analogously in reverse, hence we obtain the second part of the eigenvalue equation,

$$\begin{equation}{\mathbb{U}}_{\mathrm{o}}{\mathbb{U}}_{\mathrm{e}}\mathbf{p}={\mathbb{U}}_{\mathrm{o}}{\mathbf{p}}^{\prime }={{\Lambda}}_{\mathrm{L}}{{\Lambda}}_{\mathrm{R}}\mathbf{p}.\end{equation} \tag{ 64 }$$

Note that for z = 1, c₊ = c₋ = 0, one recovers the NESS bulk relations (40) and (44). The boundary equations given by (59)–(62) can be solved to obtain the parameters $\xi ,\omega ,z,{c}_{+},{c}_{-}\in \mathbb{C}$, together with the eigenvectors and eigenvalues of the first (leading decay-mode) orbital. Explicitly, from the left boundary equations we obtain

$$\begin{equation}\xi =\frac{z(\alpha +\beta -1)-\beta {{\Lambda}}_{\mathrm{L}}}{(\beta -1){{\Lambda}}_{\mathrm{L}}^{2}},\hspace{25.0pt}\omega =\frac{{{\Lambda}}_{\mathrm{L}}(\alpha z-{{\Lambda}}_{\mathrm{L}})}{(\beta -1)z},\hspace{25.0pt}\frac{{c}_{-}}{{c}_{+}}=\frac{{{\Lambda}}_{\mathrm{L}}^{4}}{{z}^{4}(\alpha +\beta -1)},\end{equation} \tag{ 65 }$$

and from the right ones,

$$\begin{equation}\xi =\frac{{{\Lambda}}_{\mathrm{R}}(\gamma z-{{\Lambda}}_{\mathrm{R}})}{(\delta -1)z},\hspace{25.0pt}\omega =\frac{z(\gamma +\delta -1)-\delta {{\Lambda}}_{\mathrm{R}}}{(\delta -1){{\Lambda}}_{\mathrm{R}}^{2}},\hspace{25.0pt}\frac{{c}_{+}}{{c}_{-}}=\frac{{{\Lambda}}_{\mathrm{R}}^{4}{z}^{4m+2}}{\gamma +\delta -1},\end{equation} \tag{ 66 }$$

where $m=\frac{n}{2}-2$. Pairwise identifying expressions for ξ, ω and c₊/c₋ from the left and right boundary equations, we obtain a triple of coupled equations for the unknowns Λ_L, Λ_R and z. For the details of the solution procedure and the form of the boundary vectors we refer the reader to [23]. Using this solution it can be shown that the eigenvalue of the leading decay mode Λ₁ scales with large n as ${{\Lambda}}_{1}=1-{{\Lambda}}_{1}^{\prime }/n+\mathcal{O}(1/{n}^{2})$, with

$$\begin{equation}{{\Lambda}}_{1}^{\prime }=\frac{[1-(\alpha +\beta -1)(\gamma +\delta -1)]\mathrm{log}[(\alpha +\beta -1)(\gamma +\delta -1)]}{2(\alpha +\beta -1)(\gamma +\delta -1)+\alpha \gamma -\beta \delta -2}.\end{equation} \tag{ 67 }$$

The 1/n scaling of the gap implies that the relaxation time is proportional to n, which is consistent with the ballistic transport observed in [21].

Based on these results one can formulate a conjecture for the higher orbital eigenvalues [23, 24]. Let p denote the orbital number, 1 ⩽ p ⩽ m = n/2 − 2. Then the complete set of eigenvalues Λ = Λ_LΛ_R satisfy the following set of Bethe-like equations,

$$\begin{equation}\begin{aligned}\frac{z(\alpha +\beta -1)-\beta {{\Lambda}}_{\mathrm{L}}}{(\beta -1){{\Lambda}}_{\mathrm{L}}^{2}}& =\frac{{{\Lambda}}_{\mathrm{R}}(\gamma z-{{\Lambda}}_{\mathrm{R}})}{(\delta -1)z},\\ \frac{{{\Lambda}}_{\mathrm{L}}(\alpha z-{{\Lambda}}_{\mathrm{L}})}{(\beta -1)z}& =\frac{z(\gamma +\delta -1)-\delta {{\Lambda}}_{\mathrm{R}}}{(\delta -1){{\Lambda}}_{\mathrm{R}}^{2}},\\ {(\alpha +\beta -1)}^{p}{(\gamma +\delta -1)}^{p}& ={({{\Lambda}}_{\mathrm{L}}{{\Lambda}}_{\mathrm{R}})}^{4p}{z}^{4(m-p)+2},\end{aligned}\end{equation} \tag{ 68 }$$

where p = 1 corresponds to the first orbital obtained analytically via equations (65) and (66). The eigenvalue solutions of (68) are shown in figure 6 and agree perfectly with numerical diagonalization of the Markov matrix $\mathbb{U}$. In contrast to the usual Bethe ansatz [2], this conjecture implies that we have only one momentum parameter (z) even for many-particle excitations. This is consistent with the Bethe-ansatz description of rule 54 proposed in [32], and can intuitively be understood as a consequence of a very singular dispersion relation, namely all quasi-particles (solitons) move with the same speed ±1. Another many-body effect is a large degeneracy of each eigenvalue in the orbital p which is conjectured [23] to be exponential, specifically 2^p−1.

Zoom In Zoom Out Reset image size

Expanding the last equation (68) for p = 1 in 1/m (with $m=\frac{n}{2}-2$) in the leading order of 1/n, we get that z becomes unimodular, z = e^iκ , where the momentum κ is restricted to the interval [0, π). Note that if some eigenvalues Λ_LΛ_R in the n → ∞ limit converged to points on the unit circle different from 1, this would imply persistent oscillations and/or non-decaying fluctuations in time for the model in the thermodynamic limit. Although one may naively expect that the effect of the boundaries becomes thermodynamically irrelevant, pushing all eigenvalues to collapse onto the unit circle (similarly to the closed model studied in sections 6–8), this is not the case. Rather, we have eigenvalues collapsing onto an algebraic curve (figure 6), which is inside the unit disk and only tangential to a unit circle at κ = 0 (z = 1). This indicates that the effect of the boundaries is thermodynamically nontrivial. Note that this implies that the thermodynamic limit and the long-time limit do not commute. More specifically, the stationary state is the only asymptotic state for any finite system size n, but with increasing n the time it takes to reach it diverges.

We start by introducing n₊ and n₋ to denote densities of left and right movers, given as

$$\begin{equation}{n}_{\nu }=\frac{2{N}_{\nu }}{n},\end{equation} \tag{ 92 }$$

where N_ν for ν ∈ {±} are the numbers of particles of both types on the lattice of length n (cf (33)). To account for the entropic contribution to the macrostate, we introduce ${n}_{\nu }^{\text{tot}}$ as the total density of 'slots' that can be occupied by particles of type ν ∈ {±}. In the context of thermodynamic Bethe ansatz (TBA) [4, 54] (see [55] for a pedagogical review), ${n}_{\nu }^{\text{tot}}$ corresponds to the total density (i.e. the sum of the densities of particles and holes). Classically we interpret it as the density of the effective free space in which the particles with a finite width can move [56, 57]. As a result of the interactions between particles, the total densities ${n}_{\nu }^{\text{tot}}$ exhibit a nontrivial dependence on densities n_± described by the following relation,

$$\begin{equation}{n}_{\nu }^{\text{tot}}=1-{n}_{\nu }+{n}_{-\nu },\quad \nu \in \left\{+,-\right\}.\end{equation} \tag{ 93 }$$

This identity follows from the TBA description of the quantum generalisation of the model (see [32] for the details). Intuitively, we can motivate it by realizing that two solitons of the same kind are either separated at least by two sites, or there is an oppositely moving soliton in between them, as is demonstrated by the following two diagrams,

Solitons of the same kind therefore behave as hard rods with a finite length (hence the term −n_ν ), but due to the attractive interaction with the other species (i.e. the shift after scattering is −2 sites) their total density increases and we obtain the term +n_−ν . This is consistent with the effective description in terms of the flea-gas [57], where the interaction between the particles of the same kind corresponds to the positive jump, while the scattering of oppositely-moving solitons induces the negative jump.

We restrict the discussion to the simplest class of generalised Gibbs ensembles (GGE), that are characterized by two chemical potentials, μ₊ and μ₋, corresponding to densities of right and left movers, so that the probability of a configuration $\underline{s}$ is given by

$$\begin{equation}{p}_{\underline{s}}^{\text{GGE}}=\frac{1}{{Z}_{\text{GGE}}}\enspace {\mathrm{e}}^{-{\mu }_{+}{N}_{+}(\underline{s})-{\mu }_{-}{N}_{-}(\underline{s})}.\end{equation} \tag{ 95 }$$

This GGE is equivalent to the state (29) given by the PSA with periodic boundary conditions (see section 3.3 for the details), and can be alternatively given in terms of the MPS (38) if one replaces the left and right boundary vectors with the trace (see the discussion in section 7.1).

In the large system size limit, the partition function is expressed in terms of the following phase integral involving densities of particles,

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\enspace {Z}_{\text{GGE}}=\underset{n\to \infty }{\mathrm{lim}}\int \mathrm{d}{n}_{+}\enspace \mathrm{d}{n}_{-}\enspace {\mathrm{e}}^{-\frac{n}{2}\left({\mu }_{+}{n}_{+}+{\mu }_{-}{n}_{-}-s[{n}_{+},{n}_{-}]\right)},\end{equation} \tag{ 96 }$$

where the configurational entropy density s[n₊, n₋] is given by,

$$\begin{equation}s[{n}_{+},{n}_{-}]=-\sum\limits _{\nu \in \left\{+,-\right\}}\left({n}_{\nu }\enspace \mathrm{log}\left(\frac{{n}_{\nu }}{{n}_{\nu }^{\text{tot}}}\right)+({n}_{\nu }^{\text{tot}}-{n}_{\nu })\mathrm{log}\left(1-\frac{{n}_{\nu }}{{n}_{\nu }^{\text{tot}}}\right)\right).\end{equation} \tag{ 97 }$$

In the TBA description of the quantum model, this form of entropy density is precisely the Yang–Yang entropy [32]. From the classical point of view, this form can at first sight be surprizing as one could naively expect the second term to vanish [56]. However, the particles in the model do not behave exactly as hard rods, which makes the entropy density more complicated. One can also verify the validity of (97) independently, without relying on the TBA description of the quantum model, by studying the partition sum in finite systems (see appendix C.2 for details). In the thermodynamic limit the expectation values of densities are given by evaluating this integral using the saddle-point method, which gives the following relation between n_ν and μ_ν

$$\begin{equation}\frac{{n}_{\nu }^{\text{tot}}-{n}_{\nu }}{{n}_{\nu }}={\mathrm{e}}^{{{\epsilon}}_{\nu }},\hspace{25.0pt}{{\epsilon}}_{\nu }={\mu }_{\nu }+\mathrm{log}\enspace \frac{1+{\mathrm{e}}^{-{{\epsilon}}_{\nu }}}{1+{\mathrm{e}}^{-{{\epsilon}}_{-\nu }}},\end{equation} \tag{ 98 }$$

where _ν can be interpreted as a single-particle (quasi-)energy [54]. Additionally, for later convenience, we introduce filling functions ϑ_ν as the ratios between the density and the total density,

$$\begin{equation}{\vartheta }_{\nu }=\frac{{n}_{\nu }}{{n}_{\nu }^{\text{tot}}}.\end{equation} \tag{ 99 }$$

Note that the GGE is completely determined either by the pair of chemical potentials {μ₊, μ₋}, or by the pair of filling functions {ϑ₊, ϑ₋}.

Excitations on top of these macrostates do not simply move with velocities ±1, but rather they slow down due to the interactions with other particles, and their velocities get renormalized into

$$\begin{equation}{v}_{\nu }=\nu \left(1-\frac{2{n}_{-\nu }}{1+{n}_{\nu }+{n}_{-\nu }}\right)=\frac{\nu }{1+2{\vartheta }_{-\nu }}.\end{equation} \tag{ 100 }$$

There are many equivalent ways of finding this result: one can derive it from the definition of the dressed velocity [12, 13] by using the TBA description put forward in [32]; one can obtain it by using the appropriate result for hard rods [56, 58], or an appropriate soliton gas [57, 59–61]; or alternatively the dressed velocity can be identified in terms of expectation values in the GGE [31]. Additionally, the simplicity of the result allows us to obtain it through an elementary argument introduced in [31]: let us imagine that we are tracing a right-moving soliton that starts at time t = 0 at position x = 0 and at time t it arrives to the position x. The soliton moves with the bare velocity 1 and gets displaced for two sites every time it scatters, therefore t − x = 2m_s, where m_s is the number of left-movers encountered by the tagged particle. The left-movers are moving with the dressed velocity v₋, therefore the tagged soliton will encounter all the left-movers that were at time 0 between 0 and x − v₋ t, which by identifying v₊ = x/t gives us,

$$\begin{equation}1-{v}_{+}=({v}_{+}-{v}_{-}){n}_{-},\hspace{25.0pt}1+{v}_{-}=({v}_{+}-{v}_{-}){n}_{+}.\end{equation} \tag{ 101 }$$

The second relation follows from an analogous procedure involving a tagged left-mover. Note that an additional factor of $\frac{1}{2}$ comes from the definition of left/right densities (92). This set of equations provides exactly the dressed velocities given by (100).

For any choice of the parameters {μ₊, μ₋}, the velocities v_ν are restricted to

$$\begin{equation}\left\vert {v}_{\nu }\right\vert \in \left[\frac{1}{3},1\right].\end{equation} \tag{ 102 }$$

The velocity cannot be larger than 1 (or smaller than −1), since the interactions can only slow the particles down. Alternatively, we can also explain the upper bound by the fact that the dynamics is given in terms of local mutually-commuting operators (cf (8) and (9)), which restricts the correlation spreading to the causal light-cone spreading with the speed 1 (see the upcoming discussion in section 7.4). The lower bound is at first sight less clear, but can be easily understood when we recall that two consecutive scatterings of a given soliton should be at least 3 time-steps apart, and therefore its effective speed cannot drop to 0. An example of a configuration corresponding to the soliton moving with the smallest velocity is shown in figure 9.

Zoom In Zoom Out Reset image size

Now we have all the necessary ingredients to obtain predictions for the dynamics on the Euler scale. We consider an example of the bipartitioning protocol, where at time t = 0, the left half of the system is prepared in the state determined by $({\vartheta }_{+}^{\mathrm{L}},{\vartheta }_{-}^{\mathrm{L}})$, while the state in the right half is $({\vartheta }_{+}^{\mathrm{R}},{\vartheta }_{-}^{\mathrm{R}})$. At long times, the state is locally described by a GGE parametrised by ϑ_±,ζ that slowly changes with the ray ζ = x/t,

$$\begin{equation}\underset{\begin{subarray}{c}x,t\to \infty \\ x/t=\zeta \end{subarray}}{\mathrm{lim}}\enspace {\vartheta }_{\nu }(x,t)={\vartheta }_{\nu ,\zeta }.\end{equation} \tag{ 103 }$$

The limiting values far from the junction are given by the left and right GGEs,

$$\begin{equation}\underset{\zeta \to -\infty }{\mathrm{lim}}\enspace {\vartheta }_{\nu ,\zeta }={\vartheta }_{\nu }^{\mathrm{L}},\hspace{25.0pt}\underset{\zeta \to \infty }{\mathrm{lim}}\enspace {\vartheta }_{\nu ,\zeta }={\vartheta }_{\nu }^{\mathrm{R}},\end{equation} \tag{ 104 }$$

while for the intermediate values of ζ, the filling functions satisfy the following consistency condition [12, 13],

$$\begin{equation}{\vartheta }_{\nu ,\zeta }=\begin{cases}_{\nu }^{\mathrm{L}},& {v}_{\nu }(\zeta ){ >}\zeta ,\\ {\vartheta }_{\nu }^{\mathrm{R}},& {v}_{\nu }(\zeta ){< }\zeta .\end{cases}\end{equation} \tag{ 105 }$$

Here ${v}_{\nu }(\zeta )=\frac{\nu }{1+2{\vartheta }_{-\nu ,\zeta }}$ denotes the dressed velocity in the state (ϑ_+,ζ , ϑ_−,ζ ), as given by equation (100). The consistency condition describes the profile consisting of two steps, determined by the velocity of left-movers in the state on the right and the velocity of right-movers on the left,

$$\begin{equation}({\vartheta }_{+,\zeta },{\vartheta }_{-,\zeta })=\begin{cases}({\vartheta }_{+}^{\mathrm{L}},{\vartheta }_{-}^{\mathrm{L}}),\quad & \zeta {< }-\frac{1}{1+2{\vartheta }_{+}^{\mathrm{L}}},\\ ({\vartheta }_{+}^{\mathrm{L}},{\vartheta }_{-}^{\mathrm{R}}),\quad & -\frac{1}{1+2{\vartheta }_{+}^{\mathrm{L}}}{< }\zeta {< }\frac{1}{1+2{\vartheta }_{-}^{\mathrm{R}}},\\ ({\vartheta }_{+}^{\mathrm{R}},{\vartheta }_{-}^{\mathrm{R}}),\quad & \zeta { >}\frac{1}{1+2{\vartheta }_{-}^{\mathrm{R}}}.\end{cases}\end{equation} \tag{ 106 }$$

Hydrodynamics up to the diffusive scale is described by the Navier–Stokes equation for filling functions ϑ₊, ϑ₋, [15, 62]

$$\begin{equation}{\partial }_{t}{\vartheta }_{\nu }+{v}_{\nu }{\partial }_{x}{\vartheta }_{\nu }=\sum\limits _{{\nu }_{1}}{D}_{\nu {\nu }_{1}}{\partial }_{x}^{2}{\vartheta }_{{\nu }_{1}}+\underset{\mathcal{O}\left({({\partial }_{x}\vartheta )}^{2}\right)}{\underbrace{\sum\limits _{{\nu }_{1},{\nu }_{2}}{A}_{\nu {\nu }_{1}}\left({\partial }_{x}{B}_{{\nu }_{1}{\nu }_{2}}\right){\partial }_{x}{\vartheta }_{{\nu }_{2}}}}.\end{equation} \tag{ 107 }$$

The simplicity of the TBA for RCA54 [32] allows us to find the following explicit form of the coefficients ${D}_{{\nu }_{1},{\nu }_{2}}$,

$$\begin{equation}{D}_{{\nu }_{1}{\nu }_{2}}=\frac{2{\nu }_{1}{\nu }_{2}{\vartheta }_{-{\nu }_{2}}(1-{\vartheta }_{-{\nu }_{2}})}{{(1+2{\vartheta }_{-{\nu }_{1}})}^{2}(1+2{\vartheta }_{-{\nu }_{2}})}.\end{equation} \tag{ 108 }$$

Here we obtain an additional factor of 2 w.r.t. the expressions introduced in [31, 32, 62], due to the convention used in the definition of particle densities (92). The diagonal coefficients D_ν,ν can be understood as the variance of the position of the particle of type ν due to the fluctuations in the particle density in the quasi-stationary state, and agree with the prediction of [31, 32].

For completeness, we conclude the section by reporting the coefficients ${A}_{{\nu }_{1},{\nu }_{2}}$, ${B}_{{\nu }_{1},{\nu }_{2}}$, that characterize the nonlinear term, as given by [62],

$$\begin{equation}\begin{aligned}{A}_{{\nu }_{1}{\nu }_{2}}& =\frac{(1+{\vartheta }_{+}+{\vartheta }_{-})({\delta }_{{\nu }_{1},{\nu }_{2}}-{\nu }_{1}{\nu }_{2}{\vartheta }_{{\nu }_{1}})}{1+2{\vartheta }_{-{\nu }_{1}}},\\ {B}_{{\nu }_{1}{\nu }_{2}}& ={D}_{{\nu }_{1}{\nu }_{2}}\frac{(1+2{\vartheta }_{-{\nu }_{1}})(1-{\vartheta }_{-{\nu }_{1}}+3{\vartheta }_{{\nu }_{1}})}{(1+2{\vartheta }_{+})(1+2{\vartheta }_{-})(1-{({\vartheta }_{-}+{\vartheta }_{+})}^{2})}.\end{aligned}\end{equation} \tag{ 109 }$$

Even thought the expressions are very explicit, we are not able to independently verify the validity of the quadratic term in the Navier–Stokes equation (107), since it is hard to solve it to obtain predictions in this limit.

We start by recalling that the NESS introduced in sections 3 and 4 in the absence of boundary driving (and by assuming periodic boundaries) reduces to a GGE from the previous section—see the discussion in section 3.3. We will therefore use the parametrisation of this class of stationary states in terms of the matrices W(ξ, ω), W'(ξ, ω) introduced in (42). In particular, for a periodic system of even size n,

$$\begin{equation}{\mathbf{p}}^{(n)}=\frac{1}{{Z}_{n}}\enspace \enspace \text{tr}\enspace \left({\mathbf{W}}_{1}{\mathbf{W}}_{2}^{\prime }{\mathbf{W}}_{3}\cdot \cdot \cdot {\mathbf{W}}_{n}^{\prime }\right),\hspace{25.0pt}{Z}_{n}=\enspace \text{tr}\enspace \enspace {T}^{n/2},\end{equation} \tag{ 110 }$$

where the partition function Z_n is expressed in terms of the transfer matrix T

$$\begin{equation}T=\left({W}_{0}+{W}_{1}\right)\left({W}_{0}^{\prime }+{W}_{1}^{\prime }\right).\end{equation} \tag{ 111 }$$

We recall that the chemical potentials associated to the densities of left and right movers are given by μ₊ = −log ξ and μ₋ = −log ω.

In this section we are interested in the large system-size limit, therefore we introduce the following graphical notation for the asymptotic (i.e. n → ∞) probability of finite block configurations,

where the transfer matrix ${T}^{\prime }=({W}_{0}^{\prime }+{W}_{1}^{\prime })({W}_{0}+{W}_{1})$ is obtained from T by the exchange of parameters ξ, ω. Analogously, the probabilities of configurations with odd block length can be obtained by summing over the value of the last bit,

These asymptotic probabilities can be formulated in a more explicit and convenient form by diagonalizing T and T'. However, the precise form is not essential for the upcoming discussion, and one can find it in appendix C.1.

These asymptotic probabilities exhibit a nontrivial factorization property: the conditional probability of observing kth bit to be in the state s_k , given the configuration on previous k − 1 sites (s₁, s₂, ..., s_k−1) only depends on the last three bits. Similarly, the conditional probability of s₁ given the configuration (s₂, s₃, ..., s_k ) only depends on (s₁, s₂, s₃). Explicitly, the following holds,

where we remark that this is satisfied independently of parity of the first site and the length k (i.e. it also holds when all the configurations are flipped upside-down). Note that this property, although straightforward to prove (the proof is provided in [27, 41]), is nontrivial and for example does not hold for the stationary state of a related cellular automaton [63] (rule 201; cf section 9) with otherwise similar properties.

A direct physical consequence of the factorization of asymptotic conditional probabilities is statistical independence of solitons: in the stationary state, the probability of finding a left- (or right-) moving soliton is independent of the positions of other solitons, as long as there is no left-mover (or right-mover) on the neighbouring sites. The conditional probability p_l of observing a left-mover, if we know that the neighbouring left ray does not contain a left-moving soliton is therefore a well-defined quantity and can be expressed in terms of asymptotic probabilities as follows,

where λ is the leading eigenvalue of $T=({W}_{0}+{W}_{1})({W}_{0}^{\prime }+{W}_{1}^{\prime })$ (see [27] for the details). Analogously, the conditional probability of finding a right-mover at a given site, if the neighbouring ray on the right is empty, is given by

Note that the set of probabilities (p_l, p_r) provides an equivalent parametrisation of the steady state, since the mapping (ξ, ω) → (p_l, p_r) can be inverted,

$$\begin{equation}\xi =\frac{{p}_{\mathrm{l}}(1-{p}_{\mathrm{r}})}{{(1-{p}_{\mathrm{l}})}^{2}},\hspace{25.0pt}\omega =\frac{{p}_{\mathrm{r}}(1-{p}_{\mathrm{l}})}{{(1-{p}_{\mathrm{r}})}^{2}}.\end{equation} \tag{ 117 }$$

Comparing these expressions with the TBA equations (98) and (99) we realize that p_l and p_r coincide with the filling functions ϑ₊ and ϑ₋, respectively.

We define time-states as probability distributions of time-configurations, i.e. configurations observed in time (as schematically shown in figure 11), under the assumption that the system is in a stationary state p. Explicitly, the components of the time-state q,

$$\begin{equation}{q}_{{b}_{0}\enspace {b}_{1}\enspace \cdot \cdot \cdot {b}_{m-1}}=\sum\limits _{\underline{s}\equiv ({s}_{-m},{s}_{-m+1},\cdot \cdot \cdot ,{s}_{m-1})}{p}_{\underline{s}}\enspace {\delta }_{{b}_{0},{s}_{0}^{0}}{\delta }_{{b}_{1},{s}_{1}^{1}}{\delta }_{{b}_{2},{s}_{0}^{2}}{\delta }_{{b}_{3},{s}_{1}^{3}}\cdot \cdot \cdot {\delta }_{{b}_{m-1},{s}_{\mathrm{m}\mathrm{o}\mathrm{d}(m-1,2)}^{m-1}},\end{equation} \tag{ 118 }$$

where ${s}_{k}^{t}$ is the value at site k of the configuration $\underline{s}=({s}_{-m},{s}_{-m+1},\cdot \cdot \cdot {s}_{m-1})$ propagated for t time steps started from initial data ${s}_{k}^{0}={s}_{k}$, and ${p}_{\underline{s}}$ is the asymptotic stationary probability (in our case given by equations (112) and (113)).

Zoom In Zoom Out Reset image size

To understand why for RCA54 the probability distribution q can be explicitly obtained, we imagine the following setup (schematically illustrated in figure 10): at time t = 0 we pick a random configuration, distributed according to the stationary probability distribution p, and let it evolve in time, while keeping track only of the bit at position x = 0 or x = 1 (depending on the parity of the time-step, cf figure 11). Generically, the probability of observing a soliton at time t depends on the full observed time-configuration. In our case, however, we recall that the solitons in the underlying stationary state p are statistically independent. This implies that at any time, the probability of a left (or right) mover reaching the observed site, is constant and given by p_l (or p_r), as long as there was no soliton in the previous time-step.

Explicitly, the conditional probability of observing a time-configuration (b₀, b₁, ..., b_k−1, b_k ) given the previous configuration (b₀, ..., b_k−1) therefore depends only on the last 4 bits, which implies a set of relations analogous to (114),

$$\begin{equation}\frac{{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{2k-1}}}{{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{2k-2}}}=f({b}_{2k-4},{b}_{2k-3},{b}_{2k-2},{b}_{2k-1}),\hspace{25.0pt}\frac{{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{2k}}}{{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{2k-1}}}={f}^{\prime }({b}_{2k-3},{b}_{2k-2},{b}_{2k-1},{b}_{2k}),\end{equation} \tag{ 119 }$$

where f and f' are two yet unknown functions of the last four bits. Conditional probabilities f⁽'⁾ can be straightforwardly determined by observing that short three-site time-configurations can be divided into three groups.

• (a)  
  A configuration can be inaccessible, i.e. does not appear in a time-configuration ${\underline{s}}_{x}$ corresponding to any initial configuration ${\underline{s}}^{t=0}$. To identify them, we only have to note that the following 4 configurations never appear in the local time-evolution rules (3),

  

  The probability of obtaining such configurations is 0, which means that the conditional probabilities f⁽'⁾(0, 1, 0, s) and f⁽'⁾(1, 1, 1, s) are not well defined for any s, while

  $$\begin{equation}{f}^{(\prime )}(s,0,1,0)={f}^{(\prime )}(s,1,1,1)=0.\end{equation} \tag{ 121 }$$
• (b)  
  A three-site configuration can have a unique extension into an accessible four-site configuration,

  

  These configurations represent examples, where the other choice to continue upwards would result in an inaccessible configuration. The corresponding conditional probabilities must therefore be 1,

  $$\begin{equation}{f}^{(\prime )}(s,0,1,1)={f}^{(\prime )}(0,1,1,0)=1,\quad s\in \left\{0,1\right\}.\end{equation} \tag{ 123 }$$
• (c)  
  A configuration (s₁, s₂, s₃) is valid and both (s₁, s₂, s₃, 0) and (s₁, s₂, s₃, 1) are also valid, in which case the two conditional probabilities are between 0 and 1. Explicitly, the configurations falling in this class are the following,

  

  Here the red arrows show the position in which the previously untracked (left or right moving) particle should appear for the new bit in the configuration to be equal to 1. As a consequence of statistical independence of solitons, the configurations on the bottom all arise with the conditional probability p_r or p_l, while the probabilities associated with the top configurations are 1 − p_r and 1 − p_l. Explicitly,

  $$\begin{equation}\begin{aligned}f(s,0,0,0)=f(1,1,0,0)=1-{p}_{\mathrm{r}},\quad & f(s,0,0,1)=f(1,1,0,1)={p}_{\mathrm{r}},\\ {f}^{\prime }(s,0,0,0)={f}^{\prime }(1,1,0,0)=1-{p}_{\mathrm{l}},\quad & {f}^{\prime }(s,0,0,1)={f}^{\prime }(1,1,0,1)={p}_{\mathrm{l}}.\end{aligned}\end{equation} \tag{ 125 }$$

This exhausts all the possible conditional probabilities, which allows us to express the full probability distribution q as

$$\begin{equation}{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{m-1}}={q}_{{b}_{0}{b}_{1}{b}_{2}}f({b}_{0},{b}_{1},{b}_{2},{b}_{3}){f}^{\prime }({b}_{1},{b}_{2},{b}_{3},{b}_{4})\cdot \cdot \cdot {f}^{(\prime )}({b}_{m-4},{b}_{m-3},{b}_{m-2},{b}_{m-,1}),\end{equation} \tag{ 126 }$$

where the last term in the product is f or f' for even or odd m, respectively. To completely determine q, we have to take into account the stationarity of the underlying state p, which implies the following consistency condition for any length m,

$$\begin{equation}{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{m-1}}=\sum\limits _{{b}_{-2},{b}_{-1}\in \left\{0,1\right\}}{q}_{{b}_{-2}{b}_{-1}{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{m-1}}.\end{equation} \tag{ 127 }$$

This finally provides (up to normalization) the explicit values of ${q}_{{b}_{0}{b}_{1}{b}_{2}}$ and we obtain an explicit representation of the time-state in a form similar to the patch-state ansatz introduced in section 3.2.

The time-state q can be equivalently represented in an MPA form, by encoding the appropriate factors f⁽'⁾ in the 3 × 3 matrices ${A}_{s}^{(\prime )}$,

$$\begin{equation}{A}_{0}=\left[\begin{matrix}\hfill 1-{p}_{\mathrm{r}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}{A}_{1}=\left[\begin{matrix}\hfill 0\hfill & \hfill {p}_{\mathrm{r}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}{A}_{s}^{\prime }({p}_{\mathrm{r}},{p}_{\mathrm{l}})={A}_{s}({p}_{\mathrm{l}},{p}_{\mathrm{r}}),\end{equation} \tag{ 128 }$$

and defining the appropriate boundary vectors $\left\langle L\right\vert$ and |R⟩,

$$\begin{equation}\langle L\vert =\frac{1}{1+{p}_{\mathrm{r}}+{p}_{\mathrm{l}}}\left[\begin{matrix}\hfill 1\hfill & \hfill {p}_{\mathrm{l}}\hfill & \hfill {p}_{\mathrm{r}}\hfill \end{matrix}\right],\hspace{25.0pt}\vert R\rangle =\left[\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right].\end{equation} \tag{ 129 }$$

Using these definitions, the full time-state can be succinctly summarized as,

$$\begin{equation}\mathbf{q}=\langle L\vert {\mathbf{A}}_{0}{\mathbf{A}}_{1}^{\prime }{\mathbf{A}}_{2}\cdot \cdot \cdot {\mathbf{A}}_{m-1}^{(\prime )}\vert R\rangle .\end{equation} \tag{ 130 }$$

Note that the right boundary vector |R⟩ is the right eigenvector corresponding to the eigenvalue 1 of both A₀ + A₁ and ${A}_{0}^{\prime }+{A}_{1}^{\prime }$, which immediately gives us access to the normalization of the state,

$$\begin{equation}\sum\limits _{{b}_{0},{b}_{1},\cdot \cdot \cdot ,{b}_{m-1}}{q}_{{b}_{0}{b}_{1}\cdot \cdot \cdot {b}_{m-1}}=\langle L\vert \underset{m}{\underbrace{({A}_{0}+{A}_{1})({A}_{0}^{\prime }+{A}_{1}^{\prime })\cdot \cdot \cdot ({A}_{0}^{(\prime )}+{A}_{1}^{(\prime )})}}\vert R\rangle =\langle L\vert R\rangle =1.\end{equation} \tag{ 131 }$$

The discussion in the previous subsection brings us to a related question. Let us imagine we know that for a given trajectory the time-configuration at position x is equal to ${\underline{s}}_{x}$. Is it possible to deterministically map it into the neighbouring configuration ${\underline{s}}_{x+1}$ (schematically represented in figure 12)? In other words: we wish to understand whether we can define a local and deterministic model in the space-direction, with respect to which the class of time-states introduced above is stationary. Generally, there is no guarantee that the space evolution can be expressed as a composition of local deterministic maps. In our case, however, we expect this to be possible due to the solitonic description of the model: the dynamics in space can be again understood as solitons moving in one of the two directions with a fixed velocity 1, and undergoing pairwise scattering (see e.g. figure 11). The main difference with respect to the usual dynamics in the time direction, is the nature of the interaction, which in this case temporarily speeds the particles up, rather than slowing them down.

Zoom In Zoom Out Reset image size

We therefore expect that it is possible to find a local deterministic map

$$\begin{equation}\phi :\underset{2r+1}{\underbrace{{\mathbb{Z}}_{2}{\times}{\mathbb{Z}}_{2}{\times}\cdots {\times}{\mathbb{Z}}_{2}}}\to {\mathbb{Z}}_{2},\quad r\in \mathbb{N},\end{equation} \tag{ 132 }$$

that gives the new value for a bit in the middle of the configuration of 2r + 1 sites,

$$\begin{equation}{s}_{x+1}^{t}=\phi \left({s}_{x/x-1}^{t-r},\cdot \cdot \cdot ,{s}_{x-1}^{t-2},{s}_{x}^{t-1},{s}_{x-1}^{t},{s}_{x}^{t+1},\cdot \cdot \cdot ,{s}_{x/x-1}^{t+r}\right).\end{equation} \tag{ 133 }$$

Note that χ given by (2) is of this form with r = 1. At the moment it is not clear what the support for the local space-evolution maps should be, but by examining the time-evolution rules (3), we quickly realize that the support should be larger than three (i.e. r ⩾ 2), since e.g. the last two diagrams suggest two different updated values corresponding to the time-configuration (0, 1, 1). Nonetheless, three-site time-configurations provide a starting point for construction of space-evolution maps with larger support. Again we realize that three-site configurations can be divided into three types.

• (a)  
  As we already noted earlier, time-configurations cannot contain substrings (0, 1, 0) and (1, 1, 1), as these are forbidden by the dynamic rules. The space map ϕ is therefore defined on the reduced space consisting only of allowed configurations.
• (b)  
  The first four diagrams in equation (3) already provide the deterministic space updates for three-site configurations (0, 0, 0), (1, 0, 0), (1, 0, 1) and (0, 1, 1),

  
• (c)  
  In the case of (0, 1, 1) and (1, 1, 0) one needs more information to be able to find the updated value of the central bit, since in both cases the transitions to 0 and to 1 are possible. However, it suffices to extend the configurations for two sites upwards and downwards to find a deterministic update. To demonstrate it, let us first focus on (0, 1, 1). By avoiding the forbidden configurations, there are only two ways to extend the time-configuration for two-sites at the bottom, (0, 1, 1, 0, 0) and (0, 1, 1, 0, 1),

  

  Now we are able to find a unique way to update these five-site configurations, using the rules we already have. For either one of these configurations the update rules can be applied to the bottom three sites, after which we try to set the new value at the top to either 0 and 1 and in both cases, only one of them obeys the restriction to allowed configurations,

  

  Thus we see that to find the updated value of the middle bit in the subconfiguration (0, 1, 1), it is sufficient to extend the subconfiguration two sites backwards in time. Similarly, for (1, 1, 0) we have to check the bits two sites forward in time. The rules for all such subconfigurations can be summarized by the following diagrams,

  

  where grey squares correspond to the sites that do not influence the updated value (denoted by red-bordered squares).

Together with the restriction to allowed configurations, diagrams from equations (134) and (137) completely determine the map ϕ, which has in our case support 7 (i.e. r = 3) and can be summarized as follows,

$$\begin{equation}\phi ({s}_{-3},{s}_{-2},{s}_{-1},{s}_{0},{s}_{1},{s}_{2},{s}_{3})=\begin{cases}_{1}\quad & {s}_{0}={s}_{-1}=0,\\ {s}_{-1}\quad & {s}_{0}={s}_{1}=0,\\ 0\quad & {s}_{-1}={s}_{1}=1,\\ {s}_{-3}\quad & {s}_{-1}={s}_{0}=1,\enspace {s}_{1}=0,\\ {s}_{3}\quad & {s}_{1}={s}_{0}=1,\enspace {s}_{-1}=1.\end{cases}\end{equation} \tag{ 138 }$$

Note that ϕ(s₋₃, s₋₂, s₋₁, s₀, s₁, s₂, s₃) does not depend explicitly on the values s₋₂, s₂, which means that all the maps applied in the same step (see equation (133) and figure 12) commute. The dynamics in space can be therefore understood as a composition of strictly local maps, that are in each step applied to all seven-site subconfigurations centred around the sites labelled by temporal indices of the same parity, and the order in which they are applied is not important. From this point of view, space-dynamics is, apart from the larger support of local maps, defined in perfect analogy to the usual evolution in the time-direction.

The staggered structure of time evolution, where each time-step is given in terms of mutually commuting local updates, implies a natural representation in terms of a (classical) tensor-network. The motivation for this comes from recent advances in constructing (typically nonintegrable) solvable models of many-body nonequilibrium dynamics in terms of quantum circuits, with notable examples being random [64–67] and dual-unitary circuits [68–73]. These models constitute a family of analytically tractable quantum many-body systems, for which many previously not-accessible out-of-equilibrium quantities can be obtained exactly. As we show below, RCA54 admits a convenient representation in terms of local gates, which allows us to find an equivalent diagrammatic procedure to recover the results from the previous subsections.

We start by defining two tensors represented by big and small circles, where each leg can be either in a state s = 0 or s = 1,

The big circle encodes the deterministic update (2), while the small circle forces each one of the legs to be in the same state. Using this graphical convention, the local time-evolution operator (10) can be represented as

The definition of the smaller tensors immediately implies that two circles sharing a leg can be combined into one with more legs. In particular, we note the following identity,

This allows us to combine all the local operators applied at the same time-step into one horizontal line of small and big balls positioned at sites with the opposite parity,

and ${\mathbb{U}}_{\mathrm{o}}$ takes a similar form with exchanged roles of the two tensors. Interchangeably applying ${\mathbb{U}}_{\mathrm{e}}$ and ${\mathbb{U}}_{\mathrm{o}}$ we obtain a network consisting of small and big tensors, positioned at sites with the sum of the two coordinates x + t ≡ 0 (mod 2) and x + t ≡ 1 (mod 2), respectively, as is diagrammatically shown in figure 13.

Zoom In Zoom Out Reset image size

7.4.1. Space dynamics in terms of non-deterministic dual gates

The symmetry of the tensor-network suggests a natural definition of space evolution in terms of the dual tensor $\hat{U}$, defined as,

Operators ${\hat{\mathbb{U}}}_{\mathrm{e}}$ and ${\hat{\mathbb{U}}}_{\mathrm{o}}$, describing the full one-step evolution in space, are in analogy to ${\mathbb{U}}_{\mathrm{e}}$ and ${\mathbb{U}}_{\mathrm{o}}$ given in terms of products of $\hat{U}$ centred around the sites on the dual lattice with even and odd parity respectively, as is diagrammatically shown in figure 13. One quickly realizes that $\hat{U}$ is not deterministic, which is expected, since we know that the deterministic map has the support 7. To understand how this happens, we define the following three-site projector to the space of allowed configurations,

We note that the local three-site dual operator projects to the same subspace as P,

and additionally, P and $\hat{U}$ commute if they overlap for at most one site, while they do not commute if the overlap is 2,

This allows us to show that the evolution in space can be equivalently given in terms of the operators reduced to the space of allowed configurations,

where P_e/o are defined as products of P on even/odd sites (in analogy to ${\hat{\mathbb{U}}}_{\mathrm{e}/\mathrm{o}}$). To get from the first to the second diagram, we use the relation (145), and to obtain the third diagram we note that all projectors P commute among themselves. The operators ${\tilde {\mathbb{U}}}_{\mathrm{e}}$, ${\tilde {\mathbb{U}}}_{\mathrm{o}}$ are precisely the space-evolution maps on the reduced subspace, and from the previous subsection we know that they can be expressed as a composition of seven-site deterministic maps. One can show this independently by expressing the seven-site deterministic gates in terms of the non-deterministic gate $\hat{U}$ squeezed between the relevant projectors, and then prove that the two descriptions are equivalent by utilizing appropriate local algebraic relations fulfiled by the gates and projectors. The full construction with the proof is reported in reference [28].

7.4.2. Tensor-network representation of multi-time correlation functions

To demonstrate the usefulness of the circuit formulation of RCA54, we return to the multi-time correlation functions at the same point and show an independent derivation of time-states. For simplicity we will only consider the maximum-entropy state, but a similar analysis can be done for a stationary state of the form introduced in subsection 7.1.

Before starting, we recall that by definition the expectation value of an observable a in a probability distribution on n sites p is a sum over all the configurations of the products of the probability and the value of the observable in the configurations,

$$\begin{equation}{\left\langle a\right\rangle }_{\mathbf{p}}=\sum\limits _{\underline{s}}a(\underline{s}){p}_{\underline{s}}={\left({\boldsymbol{\omega }}^{T}\right)}^{\otimes n}\cdot A\cdot \mathbf{p},\quad \boldsymbol{\omega }=\left[\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right],\end{equation} \tag{ 148 }$$

where we introduced the unnormalized one-site maximum-entropy state ω and A is a 2ⁿ × 2ⁿ diagonal matrix with the values $a(\underline{s})$. Note that the maximum entropy state on n sites corresponds to p = 2⁻ⁿ ω ^⊗n . This allows us to introduce the multi-time correlation function of one-site observables a₁, a₂, ..., a_k at times t = 0, 1, ...k − 1 as,

where represents a one-site observable and = ω. For simplicity we will assume that k ≡ 1 (mod 2), but the same can be repeated for even k. Note that these correlation functions can be directly related to time-states q as

$$\begin{equation}{C}_{{a}_{1},{a}_{2},\cdot \cdot \cdot ,{a}_{2m}}=\sum\limits _{{s}_{1},{s}_{2},\cdot \cdot \cdot ,{s}_{2m}}{q}_{{s}_{1},\cdot \cdot \cdot {s}_{2m}}\prod\limits _{j=1}^{2m}{a}_{j}({s}_{j}).\end{equation} \tag{ 150 }$$

To reduce the diagrammatic representation of ${C}_{{a}_{1},\cdot \cdot \cdot ,{a}_{k}}$, we first note that the local time-evolution operator U is deterministic and as a consequence the maximum-entropy state is invariant under it (as well as under its transpose),

Therefore the circuit (149) can be simplified by removing gates connected to ω , which reduces it to a tilted square shape. Additionally, since the observables are diagonal, they commute with the small tensor,

and we can recast the circuit in terms of dual gates $\hat{U}$ as

So far we have not taken into account any special property of RCA54 and the computational complexity of such an object would generally still grow exponentially. However, we know that in our case the complexity is constant, therefore there must be some simplification occurring.

If the gates $\hat{U}$ were deterministic, a relation analogous to (151) would be satisfied in the space direction and the correlation functions would trivialize. This does not hold precisely in our case, but the fact that the dual evolution is deterministic on the reduced subspace suggests that the gates $\hat{U}$ have some additional structure. In particular, we find that the following two few-site algebraic relations are fulfiled,

and since ${\hat{U}}^{T}=\hat{U}$, an analogous set of left-right flipped identities holds. These relations are used, together with the projector identity (145), to remove dual gates from the diagram (153) layer by layer, starting at the edges and working our way inwards. We are left with only two layers of gates, one on each side of the column of observables,

at which point the algebraic relations (154) can no longer be applied. Note that here we implicitly used the fact that observables are diagonal and therefore commute with projectors P.

By definition (cf (150)) the fully simplified tensor network (155) gives the following diagrammatic representation of components of the time-state q corresponding to the maximum entropy state,

The diagram immediately implies that the Schmidt rank does not grow with time and can be bounded from above by 4. The equivalence between this object and the corresponding time-state (130) can be straightforwardly shown by explicitly extracting 4 × 4 matrices and boundary vectors from (156). The details of the proof are reported in reference [28], together with the generalisation of the circuit representation of the multi-time correlations to the full class of Gibbs-like stationary states p.

Furthermore, analogous treatment can be straightforwardly recast for the quantum version of the model [35–37]. The main conceptual difference between the quantum and classical case is the size of the vector space, which is in the quantum case doubled (squared in dimensionality), as we are considering density matrices rather than probability distributions. Despite that, a similar set of few-site algebraic relations can be formulated, which gives access to objects analogous to (156).

The MPA representation of q gives us information about all multi-point correlation functions. In particular, one can obtain the full probability distribution of gaps between the consecutive observed particles, which in related stochastic models shows nontrivial phenomenology [53, 74, 75]. In our case, however, this behaviour is not reproduced, since the solitons passing through the origin are uncorrelated (see [27] for details).

Another quantity of interest is the one-site dynamic density–density correlation function

$$\begin{equation}C(0,t)={\left\langle \rho (t)\rho \right\rangle }_{\mathbf{p}}-{\left\langle \rho (t)\right\rangle }_{\mathbf{p}}{\left\langle \rho \right\rangle }_{\mathbf{p}},\end{equation} \tag{ 157 }$$

where ρ denotes the local density of full sites at position x = 0. The correlation function can be easily expressed in terms of the MPS (130) as,

$$\begin{equation}C(0,t)=\left\langle L\left\vert {A}_{1}{\left(({A}_{0}^{\prime }+{A}_{1}^{\prime })({A}_{0}+{A}_{1})\right)}^{t/2-1}({A}_{0}^{\prime }+{A}_{1}^{\prime }){A}_{1}\right\vert R\right\rangle -{\left\langle L\left\vert {A}_{1}\right\vert R\right\rangle }^{2},\end{equation} \tag{ 158 }$$

where we for simplicity assume even t. Since the matrices ${A}_{s}^{(\prime )}$ are finite-dimensional, the above matrix product can be easily evaluated and shown to be equal to

$$\begin{equation}C(0,t)=\frac{\left({{\Lambda}}_{1}^{\frac{t}{2}}-{{\Lambda}}_{2}^{\frac{t}{2}}\right)({({p}_{\mathrm{l}}+{p}_{\mathrm{r}})}^{2}-{p}_{\mathrm{l}}{p}_{\mathrm{r}}(1+{p}_{\mathrm{l}}+{p}_{\mathrm{r}}))+\left({{\Lambda}}_{1}^{\frac{t-2}{2}}-{{\Lambda}}_{2}^{\frac{t-2}{2}}\right){p}_{\mathrm{l}}{p}_{\mathrm{r}}({p}_{\mathrm{l}}+{p}_{\mathrm{r}})}{{(1+{p}_{\mathrm{l}}+{p}_{\mathrm{r}})}^{2}({{\Lambda}}_{1}-{{\Lambda}}_{2})},\end{equation} \tag{ 159 }$$

where Λ_1,2 are subleading eigenvalues of the transfer matrix $({A}_{0}^{\prime }+{A}_{1}^{\prime })({A}_{0}+{A}_{1})$,

$$\begin{equation}{{\Lambda}}_{1,2}=-\frac{1}{2}\left({p}_{\mathrm{l}}+{p}_{\mathrm{r}}-{p}_{\mathrm{l}}{p}_{\mathrm{r}}{\pm}\sqrt{{({p}_{\mathrm{l}}+{p}_{\mathrm{r}}-{p}_{\mathrm{l}}{p}_{\mathrm{r}})}^{2}-4{p}_{\mathrm{l}}{p}_{\mathrm{r}}}\right).\end{equation} \tag{ 160 }$$

The correlation function decays exponentially for all values of parameters, and furthermore, one can show [35, 36] that the correlations exponentially decay also when we move away from x = 0 as long as $\vert x/t\vert {< }\frac{1}{3}$. This is compatible with the hydrodynamic picture, according to which the correlations move with velocity larger than $\frac{1}{3}$ (cf (102)), and therefore we cannot observe power-law decay on rays with $\vert x/t\vert {< }\frac{1}{3}$.

We assume the finite chain with periodic boundaries of even length n, and the sites are labelled by integers in the interval [−n/2 + 1, n/2], so that the site 0 is at the middle of the chain. The length of half of the chain n/2 is assumed to be strictly larger than any finite time t we are considering. Here we use $\mathbb{U}(t)$ to denote time-evolution for the first t time-steps, i.e.

$$\begin{equation}\mathbb{U}(t)=\begin{cases}{({\mathbb{U}}_{\mathrm{o}}{\mathbb{U}}_{\mathrm{e}})}^{\frac{t}{2}},\quad & t\equiv 0\enspace (\mathrm{mod}\enspace 2),\\ {\mathbb{U}}_{\mathrm{e}}{({\mathbb{U}}_{\mathrm{o}}{\mathbb{U}}_{\mathrm{e}})}^{\frac{t-1}{2}},\quad & t\equiv 1\enspace (\mathrm{mod}\enspace 2),\end{cases}\end{equation} \tag{ 161 }$$

and in the first time-step, even sites get updated.

As we already noted in the previous section, classical observables can be understood as diagonal operators acting on the 2ⁿ -dimensional space of (macroscopic) states, so that the classical analogue of the Heisenberg picture time-evolution is given by

$$\begin{equation}a(t)=\mathbb{U}{(t)}^{-1}a\mathbb{U}(t),\end{equation} \tag{ 162 }$$

which in our case reduces to

$$\begin{equation}a(t+1)=\begin{cases}{\mathbb{U}}_{\mathrm{e}}a(t){\mathbb{U}}_{\mathrm{e}},\quad & t\equiv 0\enspace (\mathrm{mod}\enspace 2),\\ {\mathbb{U}}_{\mathrm{o}}a(t){\mathbb{U}}_{\mathrm{o}},\quad & t\equiv 1\enspace (\mathrm{mod}\enspace 2).\end{cases}\end{equation} \tag{ 163 }$$

In the following we will consider the simplest local observable: local density in the centre of the chain ρ = ρ₀, where ρ_x denotes the density at the position x,

$$\begin{equation}{\rho }_{x}={\mathbb{1}}^{\otimes n/2+x-1}\otimes \left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right]\otimes {\mathbb{1}}^{n/2-x}\equiv \left\vert 1\right\rangle {\left\langle 1\right\vert }_{x}.\end{equation} \tag{ 164 }$$

Here we introduced the shorthand symbol $\left\vert s\right\rangle {\left\langle s\right\vert }_{x}$, s = 0, 1, for the projector to one of the local physical states at site x. Additionally, we assume the following notation for the projector to the configuration (s_−m , s_−m+1, ..., s_m ) on the section [−m, m] of the lattice,

$$\begin{equation}\left\vert {s}_{-m}\cdot \cdot \cdot {s}_{m}\right\rangle \left\langle {s}_{-m}\cdot \cdot \cdot {s}_{m}\right\vert =\left\vert {s}_{-m}\right\rangle {\left\langle {s}_{-m}\right\vert }_{-m}\cdot \left\vert {s}_{-m+1}\right\rangle {\left\langle {s}_{-m+1}\right\vert }_{-m+1}\cdot \cdot \cdot \left\vert {s}_{m}\right\rangle {\left\langle {s}_{m}\right\vert }_{m}.\end{equation} \tag{ 165 }$$

Note that ρ_x (t) can be obtained from ρ(t) by a lattice spatial shift operator η,

$$\begin{equation}{\rho }_{x}(t)=\begin{cases}^{x}\rho (t){\eta }^{-x},\quad & x\equiv 0\enspace (\mathrm{mod}\enspace 2),\\ {\eta }^{x}\rho (t-1){\eta }^{-x},\quad & x\equiv 1\enspace (\mathrm{mod}\enspace 2),\end{cases}\quad {\eta }^{x}\left\vert s\right\rangle {\left\langle s\right\vert }_{y}{\eta }^{-x}=\left\vert s\right\rangle {\left\langle s\right\vert }_{y+x}.\end{equation} \tag{ 166 }$$

Furthermore, ρ_x (t) also gives immediate access to the time-evolved local density of empty sites ${\bar{\rho }}_{x}=\left\vert 0\right\rangle {\left\langle 0\right\vert }_{x}$, through the relation ${\bar{\rho }}_{x}(t)=\mathbb{1}-{\rho }_{x}(t)$. Therefore, assuming that we know ρ(t), we are in principle able to express time evolution of any local diagonal operator, since time evolution of observables is a homomorphism, (AB)(t) = A(t)B(t).

Time evolution is deterministic and local (see e.g. (151)), which implies that at time t the time-evolved observable acts nontrivially on the section of the chain between −t and t,

$$\begin{equation}\rho (t)=\sum\limits _{\underline{s}\in {\mathbb{Z}}_{2}^{2t+1}}{c}_{\underline{s}}(t)\left\vert \underline{s}\right\rangle \left\langle \underline{s}\right\vert \end{equation} \tag{ 167 }$$

where ${c}_{\underline{s}}(t)$ are the appropriate coefficients in the basis expansion. By the definition of time evolution (10) each one of these observables will be mapped into observables with a (two sites) larger support,

$$\begin{equation}\begin{aligned}& {\mathbb{U}}_{\mathrm{e}/\mathrm{o}}\left\vert {s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{t}\right\rangle \left\langle {s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{t}\right\vert {\mathbb{U}}_{\mathrm{e}/\mathrm{o}}\\ & \hspace{25.0pt}\quad =\sum\limits _{{s}_{-t-1},{s}_{t+1}\in \left\{0,1\right\}}\left\vert {s}_{-t-1}{s}_{-t}^{\prime }{s}_{-t+1}\cdot \cdot \cdot {s}_{t}^{\prime }{s}_{t+1}\right\rangle \left\langle {s}_{-t-1}{s}_{-t}^{\prime }{s}_{-t+1}\cdot \cdot \cdot {s}_{t}^{\prime }{s}_{t+1}\right\vert ,\quad \\ & \hspace{25.0pt}{s}_{j}^{\prime }=\chi ({s}_{j-1},{s}_{j},{s}_{j+1}),\end{aligned}\end{equation} \tag{ 168 }$$

where the choice between ${\mathbb{U}}_{\mathrm{e}}$ and ${\mathbb{U}}_{\mathrm{o}}$ depends on the parity of the time step (cf (163)). This immediately implies that there are exactly 4^t different configurations $\underline{s}$, referred to as accessible configurations, of length 2t + 1, for which ${c}_{\underline{s}}(t)=1$, while the other 4^t coefficients are zero. Furthermore, accessible configurations are precisely the ones for which the site at the origin was full at time t = 0, which means that the configuration $\underline{s}$ contains a soliton that was at time t = 0 at the site x = 0. The problem of time-evolution of local density can be therefore mapped onto the equivalent problem of identifying solitons in the configuration and backtracking their position in time.

By itself, this realization is not very deep and in principle the complexity of time-evolution could still grow exponentially with time t. However, in our case, the dynamics of solitons is very simple, which significantly reduces the complexity of the problem. In particular, for any soliton found in the configuration $\underline{s}=({s}_{-t},{s}_{-t+1},\cdot \cdot \cdot ,{s}_{t})$ it is possible to determine whether or not at time t = 0 it was at the origin by counting the number of scatterings it was involved in. An example of an accessible configuration is shown in figure 14. It contains several solitons, among which there is a left-mover (denoted by the red colour), that was at time t = 0 at the position x = 0. To show that the red soliton passed through the origin, one has to identify three right-moving solitons it scattered with, which are all contained in the part of the configuration to the right of it. However, one should be careful not to count too many right movers, as, e.g. the orange soliton never scattered with the red one, which is achieved by keeping track of the left-movers positioned to the right of the red soliton.

Zoom In Zoom Out Reset image size

To make this simple idea more precise, we start by splitting the coefficient ${c}_{\underline{s}}(t)$ into the sum of two terms ${c}_{\underline{s}}^{\mathrm{L}}(t)$ and ${c}_{\underline{s}}^{\mathrm{R}}(t)$, which give 1 when the soliton from the origin is moving in the left and the right direction respectively. As we will argue later, both contributions can be efficiently expressed in terms of products of matrices

$$\begin{equation}{c}_{\underline{s}}(t)=\;\;\underset{\begin{matrix}\hfill {c}_{\underline{s}}^{\mathrm{L}}(t)\hfill \end{matrix}}{\underbrace{ \left\langle {L}^{\mathrm{L}}(t)\left\vert {X}_{{s}_{-t}}^{\mathrm{L}}{Y}_{{s}_{-t+1}}^{\mathrm{L}}{X}_{{s}_{-t+2}}^{\mathrm{L}}\cdot \cdot \cdot {X}_{{s}_{t}}^{\mathrm{L}}\right\vert {R}^{\mathrm{L}}\right\rangle }}+\underset{\begin{matrix}\hfill {c}_{\underline{s}}^{\mathrm{R}}(t)\hfill \end{matrix}}{\underbrace{\left\langle {L}^{\mathrm{R}}\left\vert {X}_{{s}_{-t}}^{\mathrm{R}}{Y}_{{s}_{-t+1}}^{\mathrm{R}}{X}_{{s}_{-t+2}}^{\mathrm{R}}\cdot \cdot \cdot {X}_{{s}_{t}}^{\mathrm{R}}\right\vert {R}^{\mathrm{R}}(t)\right\rangle }},\end{equation} \tag{ 169 }$$

where ${X}_{s}^{\mathrm{L}/\mathrm{R}},{Y}_{s}^{\mathrm{L}/\mathrm{R}}\in \mathrm{E}\mathrm{n}\mathrm{d}(\mathcal{V})$, s ∈ {0, 1} are linear operators over the infinite dimensional auxiliary Hilbert space $\mathcal{V}$,

$$\begin{equation}\mathcal{V}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert c,w,n,a\rangle \vert c,w\in {\mathbb{N}}_{0},\enspace n\in \left\{0,1,2\right\},\enspace a\in \left\{0,1\right\}\right\},\end{equation} \tag{ 170 }$$

and |L^L(t)⟩, |R^L⟩, |R^R(t)⟩, |L^R⟩ boundary vectors from the auxiliary space. Since the boundary vectors change with time t, we refer to the expression (169) as the tMPA. To compactly express the matrices, we define the ladder operators c^+/−, w^+/− that change the value of c and w by one,

$$\begin{equation}\begin{aligned}{\mathbf{c}}^{+}=\sum\limits _{c,w,n,a}\left\vert c+1,w,n,a\right\rangle \left\langle c,w,n,a\right\vert ,\quad & {\mathbf{c}}^{-}={\mathbf{c}}^{+T},\\ {\mathbf{w}}^{+}=\sum\limits _{c,w,n,a}\left\vert c,w+1,n,a\right\rangle \left\langle c,w,n,a\right\vert ,\quad & {\mathbf{w}}^{-}={\mathbf{w}}^{+T},\end{aligned}\end{equation} \tag{ 171 }$$

and the projectors

$$\begin{equation}\begin{aligned}{\mathbf{c}}_{{c}_{1}{c}_{2}}=\sum\limits _{w,n,a}\left\vert {c}_{1},w,n,a\right\rangle \left\langle {c}_{2},w,n,a\right\vert ,\quad & {\mathbf{w}}_{{w}_{1}{w}_{2}}=\sum\limits _{c,n,a}\left\vert c,{w}_{1},n,a\right\rangle \left\langle c,{w}_{2},n,a\right\vert ,\\ {\mathbf{n}}_{{n}_{1}{n}_{2}}=\sum\limits _{c,w,a}\left\vert c,w,{n}_{1},a\right\rangle \left\langle c,w,{n}_{2},a\right\vert ,\quad & {\mathbf{a}}_{{a}_{1}{a}_{2}}=\sum\limits _{c,w,n}\left\vert c,w,n,{a}_{1}\right\rangle \left\langle c,w,n,{a}_{2}\right\vert .\end{aligned}\end{equation} \tag{ 172 }$$

We will elaborate on the physical interpretation of the auxiliary space in the next subsection. The operators ${X}_{s}^{\mathrm{L}}$, ${Y}_{s}^{\mathrm{L}}$ can be expressed as 3 × 3 matrices acting on the space spanned by ${\left\{\vert n\rangle \right\}}_{n=0}^{2}$ as

$$\begin{equation}\begin{aligned}{X}_{0}^{\mathrm{L}} & = \left[\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathbf{a}}_{00}+{\mathbf{a}}_{01}{\mathbf{c}}^{+}+{\mathbf{a}}_{11}{\mathbf{c}}^{+}{\mathbf{w}}^{+}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right]+{\mathbf{a}}_{11}{\mathbf{w}}_{00}\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}\\ {X}_{1}^{\mathrm{L}} & = \left[\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{a}}_{00}+{\mathbf{a}}_{01}+{\mathbf{a}}_{11}{\mathbf{w}}^{+}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{a}}_{00}+{\mathbf{a}}_{01}+{\mathbf{a}}_{11}{\mathbf{w}}^{+}\hfill \end{matrix}\right]+{\mathbf{a}}_{11}{\mathbf{w}}_{00}\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right],\\ {Y}_{0}^{\mathrm{L}} & = \left[\begin{matrix}\hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{w}}^{+}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{01}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{c}}^{+}{\left.{\mathbf{w}}^{+}\right.}^{2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{w}}^{+}\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right]+{\mathbf{a}}_{11}\left[\begin{matrix}\hfill {\mathbf{w}}_{00}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathbf{c}}^{+}{\mathbf{w}}_{10}+{\mathbf{w}}_{00}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\mathbf{w}}_{00}\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\\ {Y}_{1}^{\mathrm{L}} & = \left[\begin{matrix}\hfill 0\hfill & \hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{w}}^{+}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{c}}^{+}{\mathbf{w}}^{+}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{a}}_{00}{\mathbf{c}}^{-}{\mathbf{w}}^{+}+{\mathbf{a}}_{11}{\mathbf{c}}^{+}{\mathbf{w}}^{+}\hfill \end{matrix}\right]+{\mathbf{a}}_{11}{\mathbf{w}}_{00}\left[\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right],\end{aligned}\end{equation} \tag{ 173 }$$

while the time-dependent boundary vectors are given by

$$\begin{equation}\langle {L}^{\mathrm{L}}(t)\vert =\langle 0,t,0,0\vert ,\hspace{25.0pt}\vert {R}^{\mathrm{L}}\rangle =\sum\limits _{n=0}^{2}\vert 0,0,n,1\rangle .\end{equation} \tag{ 174 }$$

The contribution ${c}_{\underline{s}}^{\mathrm{R}}(t)$ can be obtained from ${c}_{\underline{s}}^{\mathrm{L}}(t)$ by taking a transpose, and adding additional boundary terms (see the discussion at the end of this section for details),

$$\begin{equation}\begin{aligned}{\left.{X}_{0}^{\mathrm{R}}\;\right.}^{T}& ={X}_{0}^{\mathrm{L}}-{\mathbf{c}}_{10}\left({\mathbf{a}}_{01}{\mathbf{w}}_{11}+{\mathbf{a}}_{11}({\mathbf{w}}_{10}+{\mathbf{w}}_{21})\right)\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}{\left.{X}_{1}^{\mathrm{R}}\;\right.}^{T}={X}_{1}^{\mathrm{L}},\\ {\left.{Y}_{0}^{\mathrm{R}}\;\right.}^{T}& ={Y}_{0}^{\mathrm{L}}-{\mathbf{a}}_{11}{\mathbf{c}}_{10}({\mathbf{w}}_{10}+{\mathbf{w}}_{20})\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right],\hspace{25.0pt}\langle {L}^{\mathrm{R}}\vert =\langle {R}^{\mathrm{L}}\vert + \langle 0,1,0,1\vert + \langle 0,1,2,1\vert ,\\ {\left.{Y}_{1}^{\mathrm{R}}\;\right.}^{T}& ={Y}_{1}^{\mathrm{L}}-{\mathbf{a}}_{11}{\mathbf{c}}_{10}({\mathbf{w}}_{10}+{\mathbf{w}}_{21})\left[\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right],\hspace{25.0pt}\vert {R}^{\mathrm{R}}(t)\rangle =\vert {L}^{\mathrm{L}}(t+1)\rangle .\end{aligned}\end{equation} \tag{ 175 }$$

The structure of these 3 × 3 matrices is very similar to the representation of the cubic algebra introduced in section 4.1. Indeed, operators ${X}_{s}^{\mathrm{L}}$, ${Y}_{s}^{\mathrm{L}}$, ${\left.{X}_{s}^{\mathrm{R}}\;\right.}^{T}$, ${\left.{Y}_{s}^{\mathrm{R}}\;\right.}^{T}$ can be understood as generalisation of ${W}_{s}^{(\prime )}$ (cf (42)), where the entries are no longer scalars, but rather more general operators acting on an infinite-dimensional vector space.

Before discussing the construction of tMPAs for ${c}_{\underline{s}}^{\mathrm{L}}(t)$ and ${c}_{\underline{s}}^{\mathrm{R}}(t)$, let us remark that even though the operators ${X}_{s}^{\mathrm{L}/\mathrm{R}}$, ${Y}_{s}^{\mathrm{L}/\mathrm{R}}$ are acting on an infinite-dimensional Hilbert space, we can for every finite t replace them by finite matrices, with the dimension that scales as $\mathcal{O}({t}^{2})$. In particular, the possible values of c and w in definitions of ladder operators (171) and projectors (172) can be restricted to 0 ⩽ c, w ⩽ t + 1. Note that this is an exceptional property, since typically one expects the complexity to grow exponentially with time.

8.1.1. The contribution of left-movers

We start the construction of the matrix-product representation (169) by discussing the first term ${c}_{\underline{s}}^{\mathrm{L}}(t)$. Afterwards we will show how to adapt the result to obtain also the contribution from the right-movers.

Let us consider a configuration $({s}_{-t},{s}_{-t+1},\cdot \cdot \cdot ,{s}_{t})=\underline{s}\in {\mathbb{Z}}_{2}^{2t+1}$. To figure out whether or not ${c}_{\underline{s}}^{\mathrm{L}}(t)=1$, we imagine starting at the left-most site of the configuration and moving towards the right, reading the configuration $\underline{s}$ site by site. Whenever we encounter a left-mover, we designate it the test soliton and, aiming to determine whether or not it originated from the centre, we start counting the solitons on its right. To encode the soliton counting procedure, we introduce four auxiliary degrees of freedom, |c, w, n, a⟩.

• (a)  
  The activation bit, a ∈ {0, 1}, tells us whether we are on the left (a = 0) or the right (a = 1) side of the test soliton. If the activation bit is turned off, the state splits into two parts whenever we encounter a left mover. The first part corresponds to a = 0, describing the situation in which the left mover is not the test soliton, while the second part represents the opposite case, with a = 1. In any other situation the activation bit remains unchanged.
• (b)  
  The collision counter, c ⩾ 0, represents the number of scatterings the test soliton had to undergo if it passed through the origin. As long as a = 0, the collision counter increases by 1 every two sites. If a = 1, the collision counter decreases by 1 whenever a right-mover that scattered with the test soliton is encountered. When we reach the right edge of the configuration, the value c = 0 tells us that the test soliton passed through the origin, and c ≠ 0 implies the opposite.
• (c)  
  The scattering width, w ⩾ 0, keeps track of the number of scatterings of the right-movers encountered after the test soliton. At the left edge the width w is equal to the time-step t, and every two sites it decreases by 1. Additionally, the width changes as w → w − 1 whenever a left mover on the right side of the test soliton is encountered. All the right-moving solitons that we meet after w drops to 0 could not have scattered with the test soliton.
• (d)  
  The occupation counter, n ∈ {0, 1, 2}, provides additional information about the particle content needed to appropriately change w and c. Explicitly, n = 0 if the current site is empty, n = 1 if the site is full and the left neighbour is empty, and n = 2 if the site and the left neighbour are both occupied.

At the left edge, the collision counter should be 0, and the scattering width is set to be equal to the time-step t. Furthermore, we start with the activation bit equal to 0 (since we have not yet met any test soliton), and the particle content outside of the light-cone does not influence the soliton counting, therefore we can without loss of generality assume n = 0. This implies the following form of the left boundary vector,

$$\begin{equation}\langle {L}^{\mathrm{L}}(t)\vert =\langle 0,t,0,0\vert .\end{equation} \tag{ 176 }$$

The right boundary vector should have nonzero overlap with vectors that correspond to a test soliton that passed through the origin, which is given by c = w = 0 and a = 1, while n can be arbitrary,

$$\begin{equation}\vert {R}^{\mathrm{L}}\rangle =\sum\limits _{n=0}^{2}\vert 0,0,n,1\rangle .\end{equation} \tag{ 177 }$$

To construct the matrices that correspond to the procedure summarized above, we consider 3 regimes.

Left side of the test soliton. Before the test soliton is encountered, the scattering width decreases by 1 every two sites and at the same time the collision counter should increase, while n should be appropriately adjusted to keep track of the consecutive full sites. The left action of the matrices on the sector a = 0 is therefore,

$$\begin{equation}\begin{aligned}\langle c,w,n,0\vert {X}_{s}^{\mathrm{L}}{\mathbf{a}}_{00}& =\langle c,w,s\cdot \mathrm{min}\left\{n+1,2\right\},0\vert ,\\ \langle c,w,n,0\vert {Y}_{s}^{\mathrm{L}}{\mathbf{a}}_{00}& =\langle c+1,w-1,s\cdot \mathrm{min}\left\{n+1,2\right\},0\vert .\end{aligned}\end{equation} \tag{ 178 }$$

Encountering the test soliton. Whenever a left-mover is encountered while a = 0, an additional vector with a = 1 is created. There are 4 configurations corresponding to this situation. The first two are simpler and represent a situation where a left-mover is detected while moving uninterrupted,

which can be summarized with the following matrix elements

$$\begin{equation}\langle c,w,1,0\vert {X}_{1}^{\mathrm{L}}{\mathbf{a}}_{11}=\langle c,w,2,0\vert {X}_{1}^{\mathrm{L}}{\mathbf{a}}_{11}=\langle c,w,2,1\vert .\end{equation} \tag{ 180 }$$

The other two configurations describe the soliton encountered while it is scattering,

where the grey squares denote the path of the soliton in the previous time-steps. The corresponding matrix elements are

$$\begin{equation}\langle c,w,1,0\vert {X}_{0}^{\mathrm{L}}{\mathbf{a}}_{11}=\langle c-1,w,0,1\vert ,\hspace{25.0pt}\langle c,w,1,0\vert {Y}_{0}^{\mathrm{L}}{\mathbf{a}}_{11}=\langle c,w-1,0,1\vert ,\end{equation} \tag{ 182 }$$

where in the first case we note that one scattering already occurred (hence c − 1), while in the second case the decrease of w is just due to moving to the right. This exhausts all the subconfigurations in which the bit a can flip (from the left) from 0 to 1.

Right side of the test soliton. Let us first assume w > 0. After the test soliton has been chosen, the changes of c and w are summarized as:

• (a)  
  c → c − 1 if there is a right-mover,
• (b)  
  w → w − 1 if there is a left-mover
• (c)  
  w → w − 1 every two sites.

Explicitly, there are two possible configurations of an uninterrupted right-mover appearing,

which are described by,

$$\begin{equation}\langle c,w,1,1\vert {Y}_{1}^{\mathrm{L}}=\langle c,w,2,1\vert {Y}_{1}^{\mathrm{L}}=\langle c-1,w-1,2,1\vert ,\end{equation} \tag{ 184 }$$

while the configurations in (179) imply the following,

$$\begin{equation}\langle c,w,1,1\vert {X}_{1}^{\mathrm{L}}=\langle c,w,2,1\vert {X}_{1}^{\mathrm{L}}=\langle c,w-1,2,1\vert .\end{equation} \tag{ 185 }$$

The configurations in (181) describe the two scattering solitons, therefore the corresponding matrix elements contain the decrease of both c and w,

$$\begin{equation}\begin{aligned}\langle c,w,1,1\vert {X}_{0}^{\mathrm{L}}& =\langle c-1,w-1,0,1\vert ,\\ \langle c,w,1,1\vert {Y}_{0}^{\mathrm{L}}& =(1-{\delta }_{w,0})\langle c-1,w-2,0,1\vert +{\delta }_{w,0}\langle c-1,0,0,1\vert .\end{aligned}\end{equation} \tag{ 186 }$$

All the remaining configurations do not contain any newly detected solitons,

therefore c remains constant and w decreases only due to moving two sites to the right,

$$\begin{equation}\begin{aligned}\langle c,w,0,1\vert {X}_{s}^{\mathrm{L}}=\langle c,w,s,1\vert ,\hspace{25.0pt}& \langle c,w,2,1\vert {X}_{0}^{\mathrm{L}}=\langle c,w,0,1\vert ,\\ \langle c,w,0,1\vert {Y}_{s}^{\mathrm{L}}=\langle c,w-1,s,1\vert ,\hspace{25.0pt}& \langle c,w,2,1\vert {Y}_{0}^{\mathrm{L}}=\langle c,w-1,0,1\vert .\end{aligned}\end{equation} \tag{ 188 }$$

The right-moving solitons encountered after the width drops to 0 could not scatter with the test soliton, because they were too far. Therefore c should stop decreasing,

$$\begin{equation}\langle c,0,n,1\vert {X}_{s}^{\mathrm{L}}=\langle c,0,n,1\vert {Y}_{s}^{\mathrm{L}}=\langle c,0,s\cdot \mathrm{max}\left\{2,n+1\right\},1\vert .\end{equation} \tag{ 189 }$$

Combining the matrix elements from equations (178)–(189), we finally recover the matrices given in (173). This completes the construction of the tMPA corresponding to ${c}_{\underline{s}}^{\mathrm{L}}(t)$.

8.1.2. The contribution of right-movers

The tMPA for ${c}_{\underline{s}}^{\mathrm{R}}(t)$ can be derived in a similar fashion, by simply reversing the directions of all solitons, which is achieved by exchanging the roles of the boundary vectors and by transposing the tMPA matrices. Additionally, we have to exclude all the configurations with a central right-mover that were already captured by ${c}_{\underline{s}}^{\mathrm{L}}(t)$. These are exactly the configurations where at time t = 0 the soliton in the origin is in the process of scattering, i.e. their trajectory starts with one of the following two diagrams,

In the previous case, these would correspond to the configurations with another scattering occurring exactly when w = 0. Therefore to exclude them, we have to increase w in the right boundary vector by 1, while at the same time modify the left boundary vector to allow for w not dropping to 0,

$$\begin{equation}\vert {R}^{\mathrm{R}}(t)\rangle =\vert {L}^{\mathrm{L}}(t+1)\rangle ,\hspace{25.0pt}\langle {L}^{\mathrm{R}}\vert =\langle {R}^{\mathrm{L}}\vert +\langle 0,1,0,1\vert +\langle 0,1,2,1\vert .\end{equation} \tag{ 191 }$$

Additionally, it does not suffice to map ${\left.{X}_{s}^{\mathrm{L}}\right.}^{T}\to {X}_{s}^{\mathrm{R}}$, ${\left.{Y}_{s}^{\mathrm{L}}\right.}^{T}\to {Y}_{s}^{\mathrm{R}}$, but we have to further modify them so that the following matrix elements are 0,

$$\begin{equation}\begin{aligned}{\mathbf{a}}_{11}{X}_{0}^{\mathrm{R}}\vert 1,1,1,0\rangle ={X}_{0}^{\mathrm{R}}\vert 1,1,1,1\rangle ={X}_{0}^{\mathrm{R}}\vert 1,2,1,1\rangle =0,\\ {Y}_{0}^{\mathrm{R}}\vert 1,1,1,1\rangle ={Y}_{0}^{\mathrm{R}}\vert 1,2,1,1\rangle =0,\\ {Y}_{1}^{\mathrm{R}}\vert 1,1,1,1\rangle ={Y}_{1}^{\mathrm{R}}\vert 1,2,1,1\rangle ={Y}_{1}^{\mathrm{R}}\vert 1,1,2,1\rangle ={Y}_{1}^{\mathrm{R}}\vert 1,2,2,1\rangle =0,\end{aligned}\end{equation} \tag{ 192 }$$

which gives exactly the form summarized in (175).

The tMPA provides means to study transport in the model. In this subsection, we provide two examples of physically relevant quantities that can be exactly obtained with our solution: density profile after a bipartite quench and a dynamical density–density correlation function. We finish the subsection by suggesting a formal extension of the tMPA to richer stationary states.

8.2.1. Inhomogeneous quench

Let us consider a particular bipartitioning protocol, where at time t = 0, a half-infinite chain prepared in the maximum-entropy state is joined together with a half-infinite empty chain, so that the initial state is given by,

$$\begin{equation}\mathbf{p}=\frac{1}{{2}^{n/2}}{\left[\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right]}^{\otimes n/2} \otimes {\left[\begin{matrix}\hfill 1\hfill \\ \hfill 0\hfill \end{matrix}\right]}^{\otimes n/2}.\end{equation} \tag{ 193 }$$

We wish to obtain the time-dependent profile of the particle density $\tilde {\rho }(x,t)$, which is given by the expectation value of ρ_x in the time-evolved inhomogeneous state,

$$\begin{equation}\tilde {\rho }(x,t)={\left\langle {\rho }_{x}\right\rangle }_{\mathbf{p}(t)}={\left\langle {\rho }_{x}(-t)\right\rangle }_{\mathbf{p}}=\begin{cases}_{{\eta }^{-x}\mathbf{p}},\quad & x+t\equiv 1,\\ {\left\langle \rho (t-1)\right\rangle }_{{\eta }^{-x}\mathbf{p}},\quad & x+t\equiv 0.\end{cases}\end{equation} \tag{ 194 }$$

Note that the negative sign in ρ_x (−t) together with the staggering forces a different convention with respect to (166).

The locality of time-evolution immediately implies that outside of the light-cone the density profile matches the appropriate boundary values, i.e. $\tilde {\rho }(x{< }-t,t)=\frac{1}{2}$ and $\tilde {\rho }(x{ >}t,t)=0$. For intermediate values of x, the expectation value is by definition (see equation (148) and the discussion around it)

$$\begin{align}r(\frac{t-x+1}{2},t) = {\left\langle \rho (t)\right\rangle }_{{\eta }^{-x}\mathbf{p}} & \;\;=\;\; \sum\limits _{\underline{s}}{c}_{\underline{s}}(t)\enspace {2}^{-(t-x+1)}\left[{\boldsymbol{\omega }}^{T\enspace \otimes 2t+1}\cdot \left\vert \underline{s}\right\rangle \left\langle \cdot \right\vert {\left[\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right]}^{\otimes t-x+1}\otimes {\left[\begin{matrix}\hfill 1\hfill \\ \hfill 0\hfill \end{matrix}\right]}^{\otimes t+x}\right]\\ & \enspace =\sum\limits _{{s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{-x}}{2}^{-(t-x+1)}{c}_{{s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{-x}0\cdot \cdot \cdot 0}(t),\end{align} \tag{ 195 }$$

where ω = [1 1]^T denotes a (unnormalised) one-site maximum-entropy state. The expectation value is (assuming odd x + t, cf (194)) conveniently expressed in terms of the tMPA as

$$\begin{equation}\begin{aligned}r(m,t)={2}^{-2m}& \left(\left\langle {L}^{\mathrm{L}}(t)\left\vert {\left(({X}_{0}^{\mathrm{L}}+{X}_{1}^{\mathrm{L}})({Y}_{0}^{\mathrm{L}}+{Y}_{1}^{\mathrm{L}})\right)}^{m}{X}_{0}^{\mathrm{L}}{\left({Y}_{0}^{\mathrm{L}}{X}_{0}^{\mathrm{L}}\right)}^{t-m}\right\vert {R}^{\mathrm{L}}\right\rangle \right.\\ & \left.+\left\langle {L}^{\mathrm{R}}\left\vert {\left(({X}_{0}^{\mathrm{R}}+{X}_{1}^{\mathrm{R}})({Y}_{0}^{\mathrm{R}}+{Y}_{1}^{\mathrm{R}})\right)}^{m}{X}_{0}^{\mathrm{R}}{\left({Y}_{0}^{\mathrm{R}}{X}_{0}^{\mathrm{R}}\;\right)}^{t-m}\right\vert {R}^{\mathrm{R}}(t)\right\rangle \right).\end{aligned}\end{equation} \tag{ 196 }$$

The matrices ${X}_{s}^{\mathrm{L}/\mathrm{R}}$, ${Y}_{s}^{\mathrm{L}/\mathrm{R}}$ are infinitely dimensional, therefore it is not immediately clear whether it is possible to evaluate this expression. However, their sparse structure enables us to find exact form of certain simple matrix elements (see [26] for the details). In particular, one can prove that r(m, t) takes the following form,

$$\begin{equation}\begin{aligned}r(m,t)& =\frac{3}{8}{\delta }_{m,t}+\frac{1}{8}({\delta }_{m,1}{\delta }_{t,1}+{\delta }_{m,2}{\delta }_{t,2})+\frac{1}{4}{\delta }_{m,3}{\delta }_{t,4}+{\theta }_{2m-t-3}\enspace {2}^{t-2m-1}\left(\genfrac{}{}{0.0pt}{}{2m-t-3}{t-m-1}\right)\\ & \quad +\frac{1}{8}\sum\limits _{y=t-m-2}^{2m-t-3}{2}^{-(2m-t-3)}\left(\genfrac{}{}{0.0pt}{}{2m-t-3}{y}\right)+\frac{3}{16}\sum\limits _{y=0}^{2m-t-3}{2}^{-y}\left(\genfrac{}{}{0.0pt}{}{y}{t-m-1}\right)\\ & \quad +\frac{1}{2}\sum\limits _{y=0}^{t-m-1}{2}^{-(m-1-y)}\left(\genfrac{}{}{0.0pt}{}{m-1-y}{y}\right),\end{aligned}\end{equation} \tag{ 197 }$$

where θ_x denotes the discrete Heaviside function (i.e. θ_x⩾0 = 1 and θ_x<0 = 0).

If m < (2t + 1)/3, only the last term remains and r(m, t) greatly simplifies to give the following density profile for x > −(t + 1)/3,

$$\begin{equation}\tilde {\rho }(x{ >}-\frac{t+2}{3},t)=\frac{1}{3}\left(1-{\left(-\frac{1}{2}\right)}^{\lfloor \frac{t+x+1}{2}\rfloor }\right).\end{equation} \tag{ 198 }$$

This profile consists of a front moving ballistically to the right with velocity 1 and exponentially suppressed oscillations behind it, as can also be seen in figure 15. Additionally, it exhibits an intuitive physical picture: in this section of the light-cone, all the solitons are moving to the right with the dressed velocity 1, since the right half of the lattice was empty at the beginning and there are no left-movers to slow them down. To determine the left-edge of this section, we note that the right-most left-moving solitons are moving with the velocity −1/3 (cf figure 9 and the discussion in section 6), therefore the full section of the lattice between x = −t/3 and x = t is populated only by right-movers [26].

Zoom In Zoom Out Reset image size

As we increase m (or equivalently, as we move towards the left edge of the lightcone), we are no longer able to make further simplifications to the expression (197). However, one can easily find an asymptotic expression that describes $\tilde {\rho }(x,t)$ on large scale. Using the Stirling formula it is possible to show that the density profile is well approximated by a shifted error function centred around x/t = −1/2, with the width that scales as $\sqrt{t}$,

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace \tilde {\rho }\left(-\frac{t}{2}+\zeta \sqrt{t},t\right)=\frac{1}{12}(5-\mathrm{e}\mathrm{r}\mathrm{f}(2\zeta )).\end{equation} \tag{ 199 }$$

On the rays with the constant ζ = x/t, the density profile matches the Euler scale prediction (106): the profile consists of two steps, with their velocities equal to $-\frac{1}{2}$ and 1. In the central part, the hydrodynamic prediction is given by the expectation value of density in the state determined by $({\vartheta }_{+},{\vartheta }_{-})=(\frac{1}{2},0)$, and can be shown to be equal to $\frac{1}{3}$ using the MPA representation of GGEs (see appendix C). Furthermore, the diffusive broadening around the left step is compatible with ${D}_{--}=\frac{1}{16}$, which is consistent with the expression (108) in the case ${\vartheta }_{+}=\frac{1}{2}$. The right-step does not exhibit any diffusive scaling, since the term D₊₊ vanishes for ϑ₋ = 0. The oscillations around $\frac{1}{3}$ decay exponentially with the distance from the step and are therefore not detectable on the hydrodynamic scale.

8.2.2. Dynamical structure factor

We proceed to the second example: the dynamic density–density correlation function, C(x, t), evaluated in the maximum-entropy state,

$$\begin{equation}C(x,t)=\underset{{2}^{2t+1}\tilde {C}(x,t)}{\underbrace{{\left\langle {\rho }_{0}(t){\rho }_{x}\right\rangle }_{{\mathbf{p}}_{\infty }}}}-\underset{{2}^{-2}}{\underbrace{{\left\langle \rho \right\rangle }_{{\mathbf{p}}_{\infty }}^{2}}},\hspace{25.0pt}{\mathbf{p}}_{\infty }={2}^{-n}{\left[\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right]}^{\otimes n},\end{equation} \tag{ 200 }$$

where we introduced $\tilde {C}(x,t)$ to denote the rescaled expectation value of the product of densities. Analogously to r(m, t) defined before, $\tilde {C}(x,t)$ can be expressed as a simple sum of matrix products corresponding to the left and right tMPAs,

$$\begin{equation}\begin{aligned}\tilde {C}(x,t)& =\left\langle {L}^{\mathrm{L}}(t)\left\vert {\left(({X}_{0}^{\mathrm{L}}+{X}_{1}^{\mathrm{L}})({Y}_{0}^{\mathrm{L}}+{Y}_{1}^{\mathrm{L}})\right)}^{\frac{t+x}{2}}{X}_{1}^{\mathrm{L}}{\left(({Y}_{0}^{\mathrm{L}}+{Y}_{1}^{\mathrm{L}})({X}_{0}^{\mathrm{L}}+{X}_{1}^{\mathrm{L}})\right)}^{\frac{t-x}{2}}\right\vert {R}^{\mathrm{L}}\right\rangle \\ & +\left\langle {L}^{\mathrm{R}}\left\vert {\left(({X}_{0}^{\mathrm{R}}+{X}_{1}^{\mathrm{R}})({Y}_{0}^{\mathrm{R}}+{Y}_{1}^{\mathrm{R}})\right)}^{\frac{t+x}{2}}{X}_{1}^{\mathrm{R}}{\left(({Y}_{0}^{\mathrm{R}}+{Y}_{1}^{\mathrm{R}})({X}_{0}^{\mathrm{R}}+{X}_{1}^{\mathrm{R}})\right)}^{\frac{t-x}{2}}\right\vert {R}^{\mathrm{R}}(t)\right\rangle ,\end{aligned}\end{equation} \tag{ 201 }$$

where we assumed x + t ≡ 0  (mod 2). In the opposite case, we take the advantage of the identity ${\left.C(x,t)\right\vert }_{x+t\equiv 1\enspace (\mathrm{mod}\enspace 2)}=C(x,t-1)$, and thus reduce it to the previous case. It is again possible to evaluate these sums, and we obtain the following expression [26],

$$\begin{equation}C(x,t)={2}^{-t-1}\sum\limits _{m=0}^{\frac{t-\left\vert x\right\vert -2}{2}}{4}^{m}\left(2\left(\genfrac{}{}{0.0pt}{}{t-2m-3}{m}\right)-\left(\genfrac{}{}{0.0pt}{}{t-2m-2}{m}\right)\right),\end{equation} \tag{ 202 }$$

where for simplicity of notation we assume $\left(\genfrac{}{}{0.0pt}{}{n{< }0}{k}\right)=0$ for any $k,n\in \mathbb{Z}$.

The correlation profile can be again split into two qualitatively different parts. If $\left\vert x\right\vert {< }(t+1)/3$, the spatial dependence disappears (apart from the staggering) and the correlations reduce to

$$\begin{equation}{\left.C(x,t)\right\vert }_{\left\vert x\right\vert {< }\frac{t+1}{3}}={2}^{-2t-3}\left(\left(1+\frac{\mathrm{i}}{\sqrt{7}}\right){\left(-1-\mathrm{i}\sqrt{7}\right)}^{t}+\left(1-\frac{\mathrm{i}}{\sqrt{7}}\right){\left(-1+\mathrm{i}\sqrt{7}\right)}^{t}\right).\end{equation} \tag{ 203 }$$

In this regime, correlations decay exponentially with time, $\left\vert C(0,t)\right\vert \sim {2}^{-t/2}$, and the above expression exactly matches the infinite temperature result obtained from multi-time correlation functions in section 7. Furthermore, in this regime C(x, t) is homogeneous in space and exponentially suppressed also when one considers the underlying stationary states to belong to the richer class of two-species GGEs considered in the previous sections, as a consequence of the minimal velocity of excitations being equal to 1/3 (cf (102)). This can be proven by generalising the circuit approach outlined in section 7 [35, 36].

If $\left\vert x\right\vert /t{ >}1/3$, we are left with a binomial sum that cannot be further simplified, therefore it again makes sense to find the asymptotic shape of the correlation profile. As is shown in figure 16, the profile consists of two ballistically moving peaks with velocities ${\pm}\frac{1}{2}$ that spread diffusively. Indeed, one can straightforwardly show

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace C\left({\pm}\left(\frac{t}{2}+\zeta \sqrt{t}\right),t\right)=\frac{1}{16\sqrt{t\pi }}\enspace {\mathrm{e}}^{-4{\zeta }^{2}}.\end{equation} \tag{ 204 }$$

The result nicely matches the hydrodynamics: the peaks move with dressed velocities ${v}_{{\pm}}={\pm}\frac{1}{2}$ and their diffusive broadening is governed by ${D}_{\nu \nu }=\frac{1}{16}$, ν ∈ {+, −}. Both quantities are compatible with their predictions given by (100) and (108) in the maximum-entropy state.

Zoom In Zoom Out Reset image size

8.2.3. Generalisation to generalised Gibbs states

So far we used the tMPA to express objects either evaluated in the vacuum state with only empty configuration or in the simple maximum entropy state, where all configurations are equally likely. However, the fact that the matrices ${X}_{s}^{\mathrm{L}/\mathrm{R}}$, ${Y}_{s}^{\mathrm{L}/\mathrm{R}}$ have in the space spanned by {|n⟩}_n=0,1,2 similar block structure to ${W}_{s}^{(\prime )}$, allows us to simply expand the tMPA to include the right probabilities by a simple element-wise multiplication of both sets of matrices, and thus generalise the approach to the class of GGEs given by (95). To make this point more precise, let us consider the dynamical correlation function in a state p given by the pair of parameters (ξ, ω) (corresponding to μ₊ = −log ξ, μ₋ = −log ω in (95)), which is by definition expressed as

where the state-dependent coefficient ${c}_{\underline{s}}^{{\ast}}(t)$ is defined as the product of the tMPA coefficient ${c}_{\underline{s}}(t)$ and the probability of finite configuration (113). Using the MPA formulation of these probabilities outlined in appendix C.1, the state-dependent coefficient can be rewritten as

$$\begin{equation}{c}_{{s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{t}}^{{\ast}}(t)={c}_{{s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{t}}(t)\frac{\left\langle {l}^{\prime }\left\vert {W}_{{s}_{-t}}^{\prime }{W}_{{s}_{-t+1}}\cdot \cdot \cdot {W}_{{s}_{t}}^{\prime }\right\vert r\right\rangle }{{\lambda }^{t}\left\langle {l}^{\prime }\left\vert {W}_{0}^{\prime }+{W}_{1}^{\prime }\right\vert r\right\rangle }.\end{equation} \tag{ 206 }$$

Since both the matrices describing stationary probabilities and the matrices constituting the tMPA have very similar 3 × 3 block structure (cf (173), (175), and (42)), we can in analogy to (169) express ${c}_{\underline{s}}^{{\ast}}(t)$ as a sum of the left and right time-dependent matrix products, by defining the new operators ${\left.{X}_{s}^{\mathrm{L}/\mathrm{R}}\;\right.}^{{\ast}}$, ${\left.{Y}_{s}^{\mathrm{L}/\mathrm{R}}\;\right.}^{{\ast}}$ as

$$\begin{equation}\begin{aligned}{\left.{X}_{s}^{\mathrm{L}}\right.}^{{\ast}}={X}_{s}^{\mathrm{L}}\odot {W}_{s}^{\prime },\hspace{25.0pt}& {\left.{Y}_{s}^{\mathrm{L}}\right.}^{{\ast}}=\frac{1}{\lambda }{Y}_{s}^{\mathrm{L}}\odot {W}_{s},\\ {\left.{X}_{s}^{\mathrm{R}}\;\right.}^{{\ast}}={X}_{s}^{\mathrm{R}}\odot {\left.{W}_{s}\right.}^{T} ,\hspace{25.0pt}& {\left.{Y}_{s}^{\mathrm{R}}\;\right.}^{{\ast}}=\frac{1}{\lambda }{Y}_{s}^{\mathrm{R}}\odot {\left.{W}_{s}^{\prime }\;\right.}^{T} ,\end{aligned},\end{equation} \tag{ 207 }$$

and similarly,

$$\begin{equation}\begin{aligned}\langle {\left.{L}^{\mathrm{L}}\right.}^{{\ast}}(t)\vert =\langle {L}^{\mathrm{L}}(t)\vert \odot \langle {l}^{\prime }\vert ,\hspace{25.0pt}& \langle {\left.{L}^{\mathrm{R}}\;\right.}^{{\ast}}\vert =\frac{\langle {L}^{\mathrm{R}}\vert \odot \langle r\vert }{\left\langle r\left\vert {\left.{W}_{0}^{\prime }\;\right.}^{T}+{\left.{W}_{1}^{\prime }\;\right.}^{T}\right\vert {l}^{\prime }\right\rangle },\\ \vert {\left.{R}^{\mathrm{R}}\;\right.}^{{\ast}}(t)\rangle =\vert {R}^{\mathrm{R}}(t)\rangle \odot \vert {l}^{\prime }\rangle ,\hspace{25.0pt}& \vert {\left.{R}^{\mathrm{L}}\right.}^{{\ast}}\rangle =\frac{\vert {R}^{\mathrm{L}}\rangle \odot \vert r\rangle }{\left\langle {l}^{\prime }\left\vert {W}_{0}^{\prime }+{W}_{1}^{\prime }\right\vert r\right\rangle }.\end{aligned}\end{equation} \tag{ 208 }$$

Here ⊙ denotes the element-wise multiplication, i.e. ${\left({a}_{1}\enspace {a}_{2}\enspace {a}_{3}\right)}^{T}\odot {\left({b}_{1}\enspace {b}_{2}\enspace {b}_{3}\right)}^{T}={\left({a}_{1}{b}_{1}\enspace {a}_{2}{b}_{2}\enspace {a}_{3}{b}_{3}\right)}^{T}$ (defined analogously for 3 × 3 matrices). The generalised expressions (207), (208) can be understood simply as attaching a weight ξ or ω to each left or right-mover that we observe in the soliton-counting procedure. When ξ = ω = 1 we recover the original tMPA matrices and vectors (173)–(175). Note that in the case of right-movers we have to use the transpose of the GGE matrices, and additionally the role of parameters is swapped. The explicit form of the new set of operators can be found in [41].

This construction in principle allows us to obtain C(x, t) for the full class of two-parameter GGEs (space and time translation invariant states) and, by appropriately choosing parameters, one can also obtain the tMPA formulation of the time-dependent density profile after a bipartite quench from an arbitrary pair of these GGE states. However, so far no attempt has been done to evaluate these expressions and obtain explicit results. Whether or not this is feasible is still an open question.

Even though we have been discussing RCA54 in the context of classical dynamics, there is nothing preventing us from treating it in the quantum setting. Indeed, since the one-time-step evolution operators ${\mathbb{U}}_{\mathrm{e}}$ and ${\mathbb{U}}_{\mathrm{o}}$ describe deterministic dynamics, they are also a special case of unitaries, $\left({\mathbb{U}}_{\mathrm{o}}{\mathbb{U}}_{\mathrm{e}}\right){\left({\mathbb{U}}_{\mathrm{o}}{\mathbb{U}}_{\mathrm{e}}\right)}^{{\dagger}}=\mathbb{1}$. In this respect the model represents one of the simplest interacting models, in which different measures of operator spreading can be easily studied [30, 31, 33, 34]. In particular, one such quantity is operator space entanglement entropy (OSEE), which measures the complexity of simulability of quantum dynamics in the Heisenberg picture [76–79]. In the case of RCA54 the bound of the rate of OSEE growth with time can be obtained exactly, by generalising the tMPA introduced above to the time-evolution of quantum (rather than classical) observables [33].

Time-evolution of a quantum one-site observable a can be split into four contributions, each one corresponding to a basis operator $\left\vert {s}_{1}\right\rangle \left\langle {s}_{2}\right\vert$,

$$\begin{equation}a=\left[\begin{matrix}\hfill {a}_{00}\hfill & \hfill {a}_{01}\hfill \\ \hfill {a}_{10}\hfill & \hfill {a}_{11}\hfill \end{matrix}\right],\quad \begin{aligned}a(t)& =\mathbb{U}{(t)}^{{\dagger}}a\mathbb{U}(t)={a}_{00}\mathbb{U}{(t)}^{{\dagger}}\left\vert 0\right\rangle \left\langle 0\right\vert \mathbb{U}(t)+{a}_{01}\mathbb{U}{(t)}^{{\dagger}}\left\vert 0\right\rangle \left\langle 1\right\vert \mathbb{U}(t)\\ & \quad +{a}_{10}\mathbb{U}{(t)}^{{\dagger}}\left\vert 1\right\rangle \left\langle 0\right\vert \mathbb{U}(t)+{a}_{11}\mathbb{U}{(t)}^{{\dagger}}\left\vert 1\right\rangle \left\langle 1\right\vert \mathbb{U}(t).\end{aligned}\end{equation} \tag{ 209 }$$

We immediately note that the diagonal elements,

$$\begin{equation}\mathbb{U}{(t)}^{{\dagger}}\left\vert 1\right\rangle \left\langle 1\right\vert \mathbb{U}(t)=\rho (t),\hspace{25.0pt}\mathbb{U}{(t)}^{{\dagger}}\left\vert 0\right\rangle \left\langle 0\right\vert \mathbb{U}(t)=\bar{\rho }(t)=\mathbb{1}-\rho (t),\end{equation} \tag{ 210 }$$

are given by the tMPA introduced in the previous subsections, while one of the off-diagonal terms can be expressed as the transpose of the other,

$$\begin{equation}\bar{\sigma }(t)=\mathbb{U}{(t)}^{{\dagger}}\left\vert 1\right\rangle \left\langle 0\right\vert \mathbb{U}(t)={\left(\mathbb{U}{(t)}^{{\dagger}}\left\vert 0\right\rangle \left\langle 1\right\vert \mathbb{U}(t)\right)}^{{\dagger}}=\sigma {(t)}^{{\dagger}}.\end{equation} \tag{ 211 }$$

Therefore, the only remaining thing to do is to find an efficient matrix-product description of $\sigma (t)=\mathbb{U}{(t)}^{{\dagger}}\left\vert 0\right\rangle \left\langle 1\right\vert \mathbb{U}(t)$.

Analogously to the diagonal case, at time t the time-evolved operator σ(t) nontrivially acts on the sublattice of length 2t + 1 and we can introduce coefficients ${d}_{\underline{s},\underline{b}}(t)$, such that

$$\begin{equation}\sigma (t)=\sum\limits _{\underline{s},\underline{b}\in {\mathbb{Z}}_{2}^{2t+1}}{d}_{\underline{s},\underline{b}}(t)\left\vert \underline{s}\right\rangle \left\langle \underline{b}\right\vert .\end{equation} \tag{ 212 }$$

By definition, each one of these terms will be at time t + 1 mapped into four different operators with support 2(t + 1) + 1,

$$\begin{equation}\begin{aligned}& {\mathbb{U}}_{\mathrm{e}/\mathrm{o}}\left\vert {s}_{-t}{s}_{-t+1}\cdot \cdot \cdot {s}_{t}\right\rangle \left\langle {b}_{-t}{b}_{-t+1}\cdot \cdot \cdot {b}_{t}\right\vert {\mathbb{U}}_{\mathrm{e}/\mathrm{o}}\\ & \hspace{25.0pt}=\sum\limits _{\begin{subarray}{c}{s}_{-t-1},{s}_{t+1}\\ {b}_{-t-1},{b}_{t+1}\end{subarray}}\left\vert {s}_{-t-1}{s}_{-t}^{\prime }\cdot \cdot \cdot {s}_{t}^{\prime }{s}_{t+1}\right\rangle \left\langle {b}_{-t-1}{b}_{-t}^{\prime }\cdot \cdot \cdot {b}_{t}^{\prime }{b}_{t+1}\right\vert {\delta }_{{b}_{-(t+1)},{s}_{-(t+1)}}{\delta }_{{b}_{t+1},{s}_{t+1}},\end{aligned}\end{equation} \tag{ 213 }$$

where we use the shorthand notation s_x ' = χ(s_x−1, s_x , s_x+1) and b_x ' = χ(b_x−1, b_x , b_x+1). The deterministic nature of time evolution has direct implications that can be summarised in the next three points.

• (a)  
  The coefficient ${d}_{\underline{s},\underline{b}}(t)$ can be either 0 or 1.
• (b)  
  If ${c}_{\underline{s}}(t)=0$, then also ${d}_{\underline{s},\underline{b}}(t)=0$ for all $\underline{b}$.
• (c)  
  If ${c}_{\underline{s}}(t)=1$, there exists a unique configuration $\underline{b}$, for which ${d}_{\underline{s},\underline{b}}(t)=1$. In this case, if the configurations $\underline{s}$ and $\underline{b}$ are evolved backwards in time, we will end up in a central 1 and 0 respectively, while in all intermediate time-steps, the states at the edges of the two light-cones coincide.

An example illustrating the point (c) is shown in figure 17. In this case, changing the initial bit from 1 to 0 produces an additional right-mover, therefore $\underline{s}$ can be mapped into $\underline{b}$ by adding the additional soliton and displacing all the quasi-particles that were affected by it.

Zoom In Zoom Out Reset image size

More generally, it can be shown that changing 1 to 0 in the initial configuration can be always understood as adding or removing 1 or 2 solitons, and in all these cases, there exists a simple deterministic map M_t from $\underline{s}$ to $\underline{b}$ that accounts for the relevant displacements of quasi-particles. Furthermore, as was demonstrated in reference [33], the map M_t can be expressed in terms of products of finite-dimensional matrices, with the size that does not change with time t. This immediately implies that the growth of the complexity of time-evolution is given by the size of operators constituting ${c}_{\underline{s}}(t)$, which is $\mathcal{O}({t}^{2})$. The growth of OSEE is therefore bounded from above by 2 log(t) (up to a constant factor). Note that this bound need not be saturated, and indeed, there are recent indications of OSEE increasing as $\frac{1}{2}\enspace \mathrm{log}t$ [34].

We first note that the transfer matrices T and T',

$$\begin{equation}\begin{aligned}T& =({W}_{0}+{W}_{1})({W}_{0}^{\prime }+{W}_{1}^{\prime })=\left[\begin{matrix}\hfill 1+\xi \omega \hfill & \hfill \omega \hfill & \hfill \xi \hfill \\ \hfill 1+\xi \hfill & \hfill \xi \omega \hfill & \hfill \xi \hfill \\ \hfill 1+\omega \hfill & \hfill \omega \hfill & \hfill \xi \omega \hfill \end{matrix}\right],\\ {T}^{\prime }& =({W}_{0}^{\prime }+{W}_{1}^{\prime })({W}_{0}+{W}_{1})={\left.T\right\vert }_{\xi {\leftrightarrow}\omega },\end{aligned}\end{equation} \tag{ C.1 }$$

have an isolated leading eigenvalue λ, which corresponds to the largest solution of the following cubic equation

$$\begin{equation}{\lambda }^{3}-(1+3\xi \omega ){\lambda }^{2}-(\xi +\omega +\xi \omega (1-3\xi \omega ))\lambda -\xi \omega {(1-\xi \omega )}^{2}=0.\end{equation} \tag{ C.2 }$$

The corresponding left and right leading eigenvectors can be expressed as,

$$\begin{equation}\vert l\rangle =\left[\begin{matrix}\hfill \omega \left({(\lambda -\xi \omega )}^{2}-\xi \omega \right)\hfill \\ \hfill {\omega }^{2}\left(\lambda -\xi \omega +\xi \right)\hfill \\ \hfill \xi \omega \left(\lambda -\xi \omega +\omega \right)\hfill \end{matrix}\right],\quad \vert r\rangle =\left[\begin{matrix}\hfill \omega \left(\lambda -\xi \omega +\xi \right)\hfill \\ \hfill {(\lambda -\xi \omega )}^{2}-\lambda -\xi \hfill \\ \hfill \omega \left(\lambda -\xi \omega +\omega \right)\hfill \end{matrix}\right],\quad \begin{aligned}\vert {l}^{\prime }\rangle ={\left.\vert l\rangle \right\vert }_{\xi {\leftrightarrow}\omega },\\ \vert {r}^{\prime }\rangle ={\left.\vert r\rangle \right\vert }_{\xi {\leftrightarrow}\omega },\end{aligned}\end{equation} \tag{ C.3 }$$

where ⟨l⁽'⁾|T⁽'⁾ = λ⟨l⁽'⁾|, T⁽'⁾|r⁽'⁾⟩ = λ|r⁽'⁾⟩. Since the leading eigenvalue λ is isolated, the thermodynamic expectation values (112) can be conveniently expressed as finite products with the leading eigenvectors,

The probabilities of odd configurations (113) can be expressed analogously,

where we use the fact that auxiliary space vectors (W₀ + W₁)|r⟩ and |r'⟩ are linearly dependent, i.e. there exists an $\alpha \in \mathbb{R}$ so that (W₀ + W₁)|r⟩ = α|r'⟩ (and a similar relation holds for swapped parameters).

The MPA can be used to obtain the configurational entropy of the state with fixed numbers of left and right-movers. We start by noting that the partition sum Z_n is expressed in terms of powers of the transfer matrix $T=({W}_{0}+{W}_{1})({W}_{0}^{\prime }+{W}_{1}^{\prime })$, Z_n = tr T^n/2, which can be obtained recursively by introducing the following parametrisation of T^k ,

$$\begin{equation}{T}^{k}=\left[\begin{matrix}\hfill {a}_{k}\hfill & \hfill \omega {b}_{k}\hfill & \hfill \xi {c}_{k}\hfill \\ \hfill {b}_{k}+\xi {c}_{k}\hfill & \hfill {a}_{k}-{b}_{k}\hfill & \hfill \xi {b}_{k}\hfill \\ \hfill \omega {b}_{k}+{c}_{k}\hfill & \hfill \omega {c}_{k}\hfill & \hfill {a}_{k}-{c}_{k}\hfill \end{matrix}\right].\end{equation} \tag{ C.6 }$$

By requiring T^k ⋅ T = T^k+1 we obtain a set of recursive relations that a_k , b_k and c_k have to satisfy,

$$\begin{equation}\begin{aligned}{a}_{k+1}& =(1+\xi \omega ){a}_{k}+(1+\xi )\omega {b}_{k}+(1+\omega )\xi {c}_{k},\\ {b}_{k+1}& ={a}_{k}+\xi \omega {b}_{k}+\xi {c}_{k},\\ {c}_{k+1}& ={a}_{k}+\omega {b}_{k}+\xi \omega {c}_{k}.\end{aligned}\end{equation} \tag{ C.7 }$$

From here we straightforwardly realize that Z_n can be expressed only in terms of a_n/2,

$$\begin{equation}{Z}_{n}=\enspace \text{tr}\enspace \enspace {T}^{n/2}=3{a}_{n/2}-({b}_{n/2}+{c}_{n/2})=2{a}_{n/2}+(\xi \omega -1){a}_{n/2-1},\end{equation} \tag{ C.8 }$$

and therefore Z_2k has to satisfy the same recurrence relation in k as a_k , which is obtained by reducing the set of identities of order one into a single higher-order relation,

$$\begin{equation}{Z}_{2k+2}=(1+3\xi \omega ){Z}_{2k}+\left(\xi +\omega +\xi \omega (1-3\xi \omega )\right){Z}_{2(k-1)}+\xi \omega {(\xi \omega -1)}^{2}{Z}_{2(k-2)}.\end{equation} \tag{ C.9 }$$

Now we express Z_n in terms of powers of ξ and ω,

$$\begin{equation}{Z}_{n}=\sum\limits _{x,y}{\alpha }_{n/2}(x,y){\xi }^{x}{\omega }^{y},\end{equation} \tag{ C.10 }$$

and recast the requirement (C.9) in terms of identities that have to be satisfied for coefficients α_n/2(x, y). It can be then verified that the solution takes the form

$$\begin{equation}{\alpha }_{k}(x,y)=\frac{k(k+x+y)}{(k+x-y)(k-x+y)}\left(\genfrac{}{}{0.0pt}{}{k+x-y}{y}\right)\left(\genfrac{}{}{0.0pt}{}{k-x+y}{x}\right),\end{equation} \tag{ C.11 }$$

where we implicitly assume α_k (x, y) = 0 for k ⩽ ±(x − y). Asymptotically, the combinatorial term reduces to the configurational entropy introduced in (97),

$$\begin{equation}\underset{n\to \infty }{\mathrm{lim}}\left({\alpha }_{n/2}\left(\frac{n}{2}{n}_{+},\frac{n}{2}{n}_{-}\right)-{\mathrm{e}}^{\frac{n}{2}s[{n}_{+},{n}_{-}]}\right)=0.\end{equation} \tag{ C.12 }$$