Exact matrix product decay modes of a boundary driven cellular automaton

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

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

The rule is given by a deterministic local mapping on a diamond-shaped plaquette $\newcommand{\ZZ}{\mathbb{Z}} \chi : \ZZ_2\times \ZZ_2\times \ZZ_2 \to \ZZ_2$. A south site s_S is determined by a north, west and east sites

$$\begin{align} \newcommand{\ZZ}{\mathbb{Z}} \displaystyle s_{\rm S} = \chi({s_{\rm W},s_{\rm N},s_{\rm E}}) = s_{\rm N} + s_{\rm W} + s_{\rm E} + s_{\rm W} s_{\rm E} \pmod{2}, \quad s_{\rm S},s_{\rm N},s_{\rm W},s_{\rm E} \in \ZZ_2.\label{chi} \nonumber \end{align} \tag{ 1 }$$

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

We thus proceed to formulate the time evolution of the full probability distribution $p_{(s_1, s_2, \ldots, s_n)}$, which we call a state³. The state space is a convex subset of a vector space $\newcommand{\ot}[1]{{\widetilde{#1}}} \newcommand{\RR}{\mathbb{R}} \newcommand{\R}{{\rm{R}}} {{\mathcal S}} = \RR^{{\mathcal C}} = (\RR^2){\hspace{0pt}}^{\otimes n}$ of probability distributions over configurations ${\bf p}=(\,p_0, p_1, \ldots, p_{2^n-1}) \in {{\mathcal S}}$ satisfying the non-negativity and normalization conditions, $p_s \geqslant 0$, $\sum_{s=0}^{2^n-1} p_s = 1$. Here, s is a binary coded configuration $s=\sum_{k=1}^n 2^{n-k} s_k$. The local Rule 54 can then be given in terms of a three-site $2^3 \times 2^3$ permutation matrix P

$$\begin{align} \displaystyle P_{(s,s^{\prime},s^{\prime\prime}),(t,t^{\prime},t^{\prime\prime})} = \delta_{s,t} \delta_{s^{\prime},\chi(t,t^{\prime},t^{\prime\prime})} \delta_{s^{\prime\prime},t^{\prime\prime}}, \nonumber \end{align} \tag{ 2 }$$

or

$$\begin{align*} \displaystyle P = \left(\begin{array}{@{}cccccccc@{}} \, 1\, & & & & & & & \nonumber \\ & & & \, 1\, & & & & \nonumber \\ & & \, 1\, & & & & & \nonumber \\ & \,1\, & & & & & & \nonumber \\ & & & & & & \,1\, & \nonumber \\ & & & & & & & \,1\, \nonumber \\ & & & & \,1\, & & & \nonumber \\ & & & & & \,1\, & & \end{array}\right), \nonumber \end{align*}$$

such that $\newcommand{\one}{\mathbb {1}} P^2=\one$. The local update rule is in turn embedded into ${\rm End}({{\mathcal S}})$ as $\newcommand{\ot}[1]{{\widetilde{#1}}} \newcommand{\one}{\mathbb {1}} P_{k, k+1, k+2} = \one_{2^{k-1}}\otimes P \otimes \one_{2^{n-k-2}}$ acting on any triple of neighboring sites $k, k+1, k+2$. The time evolution of the state vector ${\bf p}(t) \in {{\mathcal S}}$ starting from some initial state ${\bf p}(0)$ is written as

$$\begin{align} \displaystyle {\bf p}(t)=U^t {\bf p}(0), \nonumber \end{align} \tag{ 3 }$$

where

$$\begin{align} \displaystyle U=U_{\rm o} U_{\rm e}, \label{U} \nonumber \end{align} \tag{ 4 }$$

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

$$\begin{align} \displaystyle U_{\rm e}=P_{123}P_{345}\cdots P_{n-3,n-2,n-1} P^{\rm{R}}_{n-1,n}, \label{Ue} \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \displaystyle U_{\rm o}= P^{\rm{L}}_{12} P_{234} P_{456} \cdots P_{n-2,n-1,n}. \label{Uo} \nonumber \end{align} \tag{ 6 }$$

The boundary propagators

$$\begin{align} \newcommand{\one}{\mathbb {1}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle P^{\rm{L}}_{12}=P^{\rm{L}} \otimes \one_{2^{n-2}}, \qquad P^{\rm{R}}_{n-1,n}= \one_{2^{n-2}} \otimes P^{\rm{R}}, \nonumber \end{align} \tag{ 7 }$$

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

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

$$\begin{align} \displaystyle [P_{k-1,k,k+1},P_{k^{\prime}-1,k^{\prime},k^{\prime}+1}] = 0,\quad {\rm if}\quad |k-k^{\prime}| \geqslant 2. \nonumber \end{align} \tag{ 8 }$$

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

$$\begin{align} \displaystyle [P^{\rm L}_{12},P_{234}]=0,\qquad [P_{123},P^{\rm R}_{34}]=0, \label{boundcomm} \nonumber \end{align} \tag{ 9 }$$

so the order of factors in either U_e or U_o, equations (5) and (6), is irrelevant.

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle U {\bf p}=\Lambda {\bf p}, \nonumber \end{align} \tag{ 10 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle U_{\rm e} {\bf p}=\Lambda_{\rm{L}} {\bf p^{\prime}}, \qquad U_{\rm o} {\bf p^{\prime}}=\Lambda_{\rm{R}} {\bf p}, \qquad \Lambda=\Lambda_{\rm{R}} \Lambda_{\rm{L}}. \label{eigensplit} \nonumber \end{align} \tag{ 11 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle {\bf p}(t) = {\bf p}_0 + \sum_{j\geqslant 1} c_j \Lambda_j^t {\bf p}_j, \nonumber \end{align} \tag{ 12 }$$

where c_j are appropriate constants depending on the initial state.

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

Zoom In Zoom Out Reset image size

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

1.1. Conditional boundary driving

For the conditional driving the boundary stochastic matrices read

$$\begin{align} \displaystyle &P^{\rm{L}}=\left(\begin{array}{@{}cccc@{}} \alpha & 0 & \alpha & 0 \nonumber \\ 0 & \beta & 0 & \beta \nonumber \\ 1-\alpha & 0 & 1-\alpha & 0 \nonumber \\ 0 & 1-\beta & 0 & 1-\beta \nonumber \\ \end{array} \right), \quad P^{\rm{R}}=\left(\begin{array}{@{}cccc@{}} \gamma & \gamma & 0 & 0 \nonumber \\ 1-\gamma & 1-\gamma & 0 & 0 \nonumber \\ 0 & 0 & \delta & \delta \nonumber \\ 0 & 0 & 1-\delta & 1-\delta \nonumber \\ \end{array} \right), \label{cond} \nonumber \end{align} \tag{ 13 }$$

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

1.2. Bernoulli boundary driving

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

$$\begin{align} \displaystyle &P^{\rm{L}}=\left(\begin{array}{@{}cccc@{}} \frac{1}{2} & 0 & \frac{1}{2} & 0 \nonumber \\ 0 & 1-\alpha & 0 & \beta \nonumber \\ \frac{1}{2} & 0 & \frac{1}{2} & 0 \nonumber \\ 0 & \alpha & 0 & 1-\beta \nonumber \\ \end{array} \right), \qquad P^{\rm{R}}=\left(\begin{array}{@{}cccc@{}} \frac{1}{2} & \frac{1}{2} & 0 & 0 \nonumber \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 \nonumber \\ 0 & 0 & 1-\gamma & \delta \nonumber \\ 0 & 0 & \gamma & 1-\delta \nonumber \\ \end{array} \right), \label{bern} \nonumber \end{align} \tag{ 14 }$$

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

$$\begin{align*} \displaystyle \mu:=\alpha-\beta, \qquad \nu:=\beta-1, \qquad \sigma:=\gamma-\delta, \qquad \rho:=\delta-1. \nonumber \end{align*}$$

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

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

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

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

We begin by defining a vector

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle {\bf W}= \left(\begin{array}{@{}c@{}} W_0(\xi, \omega) \nonumber \\ W_1(\xi, \omega) \end{array}\right), \nonumber \end{align} \tag{ 15 }$$

with components W_s being $4\times 4$ matrices

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle W_0(\xi, \omega)=\left(\begin{array}{@{}cccc@{}} 1 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \xi & \xi & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \end{array} \right), \qquad W_1(\xi, \omega)=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & \xi & 1 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 1 & \omega \nonumber \\ \end{array} \right), \label{Wrepr} \nonumber \end{align} \tag{ 16 }$$

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

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle {\bf W}^{\prime}(\xi, \omega)={\bf W}(\omega,\xi). \nonumber \end{align} \tag{ 17 }$$

The key property explored in this paper is a simple three-site cubic algebraic relation⁵ which shall provide a cancellation mechanism to be used later for constructing the eigenstates of the Markov matrix

$$\begin{align} \displaystyle P_{123} {\bf W}_1 S\, {\bf W}_2 {\bf W}^{\prime}_3={\bf W}_1 {\bf W}^{\prime}_2 {\bf W}_3 S, \label{bulk1} \nonumber \end{align} \tag{ 18 }$$

where S is a constant 'delimiter' matrix,

$$\begin{align} \newcommand{\one}{\mathbb {1}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle S=\one \otimes \sigma^{\rm{x}}=\left(\begin{array}{@{}cccc@{}} 0 & 1 & 0 & 0 \nonumber \\ 1 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ \end{array} \right), \label{Sdef} \nonumber \end{align} \tag{ 19 }$$

satisfying $\newcommand{\one}{\mathbb {1}} S^2=\one_4$. By interchanging ξ and ω in (18), i.e. interchanging ${\bf W}$ and ${\bf W}^{\prime}$, and multiplying with $P_{123}$ (noting that $\newcommand{\one}{\mathbb {1}} P^2=\one_8$) we obtain a dual bulk relation

$$\begin{align} \displaystyle P_{123} {\bf W}^{\prime}_1 {\bf W}_2 {\bf W}^{\prime}_3 S = {\bf W}^{\prime}_1 S\, {\bf W}^{\prime}_2 {\bf W}_3. \label{bulk1b} \nonumber \end{align} \tag{ 20 }$$

Note that equation (18) in fact represents 8 matrix product identities, $W_s S W_{\chi(ss^{\prime}s^{\prime\prime})} W^{\prime}_{s^{\prime\prime}} = W_s W^{\prime}_{s^{\prime}} W_{s^{\prime\prime}} S$, and analogously for equation (20), by writing out physical space components $s, s^{\prime}, s^{\prime\prime}\in\{0, 1\}$ for a vector on a triple of consecutive physical sites (denoted in (18) and (20) as 123). Equations (18) and (20) thus represent a set of algebraic relations among $W_0, W_1, W^{\prime}_0, W^{\prime}_1, S$ which can be straightforwardly verified for the representation (16).

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

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p}=\bra{{\bf l}_1}\mmu{W}{2}\mmup{W}{3} \mmu{W}{4} \mmup{W}{5}\cdots\mmup{W}{n-3}\mmu{W}{n-2}\ket{{\bf r}_{n-1,n}}, \label{pness} \nonumber \end{align} \tag{ 21 }$$

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p^{\prime}}=\bra{{\bf l}^{\prime}_{12}}\mmu{W}{3}\mmup{W}4\mmu{W}{5}\mmup{W}{6}\cdots\mmup{W}{n-2}\mmu{W}{n-1}\ket{{\bf r}^{\prime}_{n}}. \label{ppness} \nonumber \end{align} \tag{ 22 }$$

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

$$\begin{align} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle P^{\rm{L}}_{12} \bra{{\bf l}^{\prime}_{12}}=\lambda_{\rm{L}} \bra{{\bf l}_1} \mmu{W}{2} S, \label{bound1a} \nonumber \end{align} \tag{ 23 }$$

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle P_{123} \bra{{\bf l}_{1}} \mmu{W}{2} \mmup{W}{3}=\bra{{\bf l}^{\prime}_{12}} \mmu{W}{3}S,\label{bound2a} \nonumber \end{align} \tag{ 24 }$$

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle P^{\rm R}_{12} \ket{{\bf r}_{12}}=\mmup{W}{1}S\ket{{\bf r}^{\prime}_2}, \label{bound3a} \nonumber \end{align} \tag{ 25 }$$

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle P_{123} \mmup{W}{1} \mmu{W}{2}\ket{{\bf r}^{\prime}_{3}} =\lambda_{\rm{R}} \mmup{W}{1}S\ket{{\bf r}_{23}}. \label{bound4a} \nonumber \end{align} \tag{ 26 }$$

Specifically, writing out $U_{\rm e} {\bf p}$ in terms of (5) and the ansatz (21), one first uses equation (25) in order to introduce the delimiter S in a string $\cdots {\bf W}^{\prime}_{n-3}{\bf W}_{n-2}{\bf W}^{\prime}_{n-1}S$ and then implements $\cdots P_{n-5, n-4, n-3}P_{n-3, n-2, n-1}$ via the dual bulk relation (20) in order to move the delimiter S across to the left end where it is then absorbed, via $\newcommand{\one}{\mathbb {1}} S^2=\one$, by applying another boundary equation (24), arriving to $U_{\rm e} {\bf p}=\lambda_{\rm R}{\bf p}^{\prime}$. Analogously we proceed with $U_{\rm o} {\bf p}^{\prime}$, in terms of (6) and ansatz (22), now implementing the boundary equations (23) and (26) to carry S from left to right via the bulk relation (18), ending with $U_{\rm e} {\bf p}^{\prime}=\lambda_{\rm L}{\bf p}$. Thus, ${\bf p}$ is an eigenvector of U with an eigenvalue⁶ $\lambda=\lambda_{\rm L}\lambda_{\rm R}$.

As for an explicit example, consider $n=6$ sites and even part of the propagator U_e

$$\begin{align} \newcommand{\mmup}[2]{{{\bf #1}^{\prime}_{#2}}} \newcommand{\mmu}[2]{{{\bf #1}_{#2}}} \newcommand{\mm}[1]{{{\bf #1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle U_{\rm e} {\bf p} &=P_{123}P_{345} P^{\rm{R}}_{56} \bra{{\bf l}_1}\mmu{W}{2}\mmup{W}{3} \mmu{W}{4} \ket{{\bf r}_{56}} \nonumber \\ &= P_{123}P_{345} \bra{{\bf l}_1}\mmu{W}{2}\mmup{W}{3} \mmu{W}{4} \mmup{W}{5}S\ket{{\bf r}^{\prime}_6} \nonumber \\ &= P_{123} \bra{{\bf l}_1}\mmu{W}{2}{\bf W}^{\prime}_3 S {\bf W}^{\prime}_4 {\bf W}_5\ket{{\bf r}^{\prime}_6} \nonumber \\ &= \bra{{\bf l}^{\prime}_{12}} \mmu{W}{3}S^2 {\bf W}^{\prime}_4 {\bf W}_5\ket{{\bf r}^{\prime}_6} \nonumber \\ &={\bf p^{\prime}}. \nonumber \end{align} \tag{ 27 }$$

Solving the full set of boundary equations (23)–(26) should fix all the unknown parameters in the MPA (21) and (22) as the bulk relations are automatically satisfied. The solution is unique up to an irrelevant transformation of boundary vectors. We first focus on the conditional driving case (13) and then state the results for Bernoulli driving (14).

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

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{(\alpha +\beta -1)-\beta \lambda_{\rm{L}}}{(\beta -1) \lambda_{L}^{2}}, \label{nessspc1} \nonumber \end{align} \tag{ 28 }$$

$$\begin{align} \displaystyle \omega =\frac{\lambda_{\rm{L}} (\alpha -\lambda_{\rm{L}})}{(\beta -1)}, \label{nessspc2} \nonumber \end{align} \tag{ 29 }$$

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

$$\begin{align} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle &\bra{l_0}=\left(0, \frac{1}{\lambda_{\rm{L}}}, \frac{\alpha (\alpha +\beta -1-\lambda_{\rm{L}})}{\beta \lambda_{\rm{L}}^2 (\alpha -\lambda_{\rm{L}})},1 \right), \nonumber \\ &\bra{l_1}= \left(0, \frac{1-\beta }{\beta \lambda_{\rm{L}}}, \frac{(1-\alpha) (\alpha +\beta-1-\lambda_{\rm{L}})}{\beta \lambda_{\rm{L}}^2 (\alpha -\lambda_{\rm{L}})},\frac{1}{\beta }-1\right), \nonumber \\ &\bra{l^{\prime}_{0,0}}=\left(\frac{\alpha (\alpha +\beta -1-\lambda_{\rm{L}})}{(\beta -1) \beta \lambda_{\rm{L}}}, \frac{\alpha +\beta -1-\lambda_{\rm{L}}}{\beta -1}, 0, 0 \right), \nonumber \\ &\bra{l^{\prime}_{0,1}}=\left(0, 0, \frac{\lambda_{\rm{L}} (\alpha -\lambda_{\rm{L}}) (\alpha +\beta-1-\lambda_{\rm{L}})}{(\beta -1) (\alpha+\beta -1 -\beta \lambda_{\rm{L}})},\frac{\alpha (\alpha +\beta -1-\lambda_{\rm{L}})}{(\beta -1) \beta \lambda_{\rm{L}}} \right), \nonumber \\ &\bra{l^{\prime}_{1,0}}=\left(\frac{(\alpha -\lambda_{\rm{L}}) (\alpha +\beta-1 -\lambda_{\rm{L}})}{(1-\beta) \beta \lambda_{\rm{L}}}, -\frac{\alpha +\beta-1 -\lambda_{\rm{L}}}{\beta }, 0, 0 \right), \nonumber \\ &\bra{l^{\prime}_{1,1}}= \left(0, 0, \frac{(1-\alpha) \lambda_{\rm{L}} (\alpha +\beta -1-\lambda_{\rm{L}})}{\beta (\alpha+\beta-1 -\beta \lambda_{\rm{L}})},\frac{(1-\alpha)(\alpha +\beta-1 -\lambda_{\rm{L}})}{(\beta -1) \beta \lambda_{\rm{L}}} \right). \nonumber \end{align} \tag{ 30 }$$

The right boundary equations (25) and (26), on the other hand, give the following unique solutions for the spectral parameters

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\lambda_{\rm{R}} (\gamma -\lambda_{\rm{R}})}{(\delta -1) }, \label{nessspc3} \nonumber \end{align} \tag{ 31 }$$

$$\begin{align} \displaystyle \omega =\frac{(\gamma +\delta -1)-\delta \lambda_{\rm{R}}}{(\delta -1) \lambda_{\rm{R}}^2},\label{nessspc4} \nonumber \end{align} \tag{ 32 }$$

and the the right boundary vectors

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle &\ket{r_{0,0}}=\left(\frac{\gamma (\gamma +\delta-1 -\lambda_{\rm{R}})}{\delta (\gamma +\delta -1 -\delta \lambda_{\rm{R}})}, 0, \frac{\gamma +\delta-1 -\lambda_{\rm{R}}}{(\delta -1) \lambda_{\rm{R}}}, 0 \right), \nonumber \\ &\ket{r_{0,1}} =\left(\frac{(\lambda_{\rm{R}}-\gamma) (\gamma +\delta-1 -\lambda_{\rm{R}})}{\delta (\gamma +\delta-1 - \delta \lambda_{\rm{R}})}, 0,- \frac{\gamma +\delta-1 -\lambda_{\rm{R}}}{\delta \lambda_{\rm{R}}}, 0 \right), \nonumber \\ &\ket{r_{1,0}}= \left(0, \frac{(\gamma -\lambda_{\rm{R}}) (\gamma +\delta -1-\lambda_{\rm{R}})}{(\delta-1) (\gamma+\delta -1 -\delta \lambda_{\rm{R}})}, 0, \frac{\gamma \lambda_{\rm{R}} (\gamma -\lambda_{\rm{R}}) (\gamma +\delta-1 -\lambda_{\rm{R}})}{(\delta -1) \delta (\gamma+\delta -1 -\delta \lambda_{\rm{R}})} \right), \nonumber \\ &\ket{r_{1,1}}=\left(0, \frac{(1-\gamma) (\gamma +\delta-1 -\lambda_{\rm{R}})}{\delta (\gamma+\delta -1 -\delta \lambda_{\rm{R}})}, 0, -\frac{(\gamma -1) \lambda_{\rm{R}} (\gamma-\lambda_{\rm{R}}) (\gamma +\delta-1 -\lambda_{\rm{R}})}{(\delta -1) \delta (\gamma+\delta -1 -\delta \lambda_{\rm{R}})} \right), \nonumber \\ &\ket{r^{\prime}_0}=\left(0, \frac{\gamma \lambda_{\rm{R}} (\gamma +\delta-1 -\lambda_{\rm{R}})}{\delta (\gamma +\delta -1-\delta \lambda_{\rm{R}})}, 1, \frac{\lambda_{\rm{R}}^2(\gamma -\lambda_{\rm{R}})}{\gamma+\delta -1 -\delta \lambda_{\rm{R}}}\right), \nonumber \\ &\ket{r^{\prime}_1}=\left(0, \frac{(1-\gamma) \lambda_{\rm{R}} (\gamma +\delta-1 -\lambda_{\rm{R}})}{\delta (\gamma +\delta -1-\delta \lambda_{\rm{R}})},\frac{1}{\delta }-1, \frac{(1-\delta) \lambda_{\rm{R}}^2 (\gamma -\lambda_{\rm{R}})}{\delta (\gamma+\delta -1-\delta \lambda_{\rm{R}})} \right). \nonumber \end{align} \tag{ 33 }$$

Now in order to get a consistent solution on both the left and right boundary we demand the two pairs of the spectral parameters ξ and ω in equations (28), (29), (31) and (32) to be equal. This gives us a closed pair of equations for $\lambda_{\rm L}$ and $\lambda_{\rm R}$,

$$\begin{align} \displaystyle \frac{(\alpha +\beta -1)-\beta \lambda_{\rm{L}}}{(\beta -1) \lambda_{L}^{2}}=\frac{\lambda_{\rm{R}} (\gamma -\lambda_{\rm{R}})}{(\delta -1) }, \nonumber \end{align} \tag{ 34 }$$

$$\begin{align} \displaystyle \frac{\lambda_{\rm{L}} (\alpha -\lambda_{\rm{L}})}{(\beta -1)}=\frac{(\gamma +\delta -1)-\delta \lambda_{\rm{R}}}{(\delta -1) \lambda_{\rm{R}}^2}. \nonumber \end{align} \tag{ 35 }$$

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

$$\begin{align} \displaystyle &(\lambda-1) \label{char1}\nonumber \\ &\times\left(\lambda ^3+\lambda ^2 (1-\alpha \gamma)+\lambda [\beta \delta -(\alpha +\beta -1) (\gamma +\delta -1)]-(\alpha +\beta -1) (\gamma +\delta -1)\right)=0. \nonumber \end{align} \tag{ 36 }$$

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

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

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\lambda_{\rm{L}}+\mu \left(1-\lambda_{\rm{L}}\right)}{(\mu +1) \lambda_{\rm{L}}^2}, \nonumber \end{align} \tag{ 37 }$$

$$\begin{align} \displaystyle \omega=\frac{\lambda_{\rm{L}} \left(2 \lambda_{\rm{L}}-1\right)}{\mu +1}, \nonumber \end{align} \tag{ 38 }$$

and on the right,

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\lambda_{\rm{R}} \left(2 \lambda _{\rm{R}}-1 \right)}{\sigma +1}, \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \displaystyle \omega=\frac{\lambda_{\rm{R}}+\sigma \left(1-\lambda_{\rm{R}}\right)}{(\sigma+1) \lambda_{\rm{R}}^2}, \nonumber \end{align} \tag{ 40 }$$

yielding the characteristic polynomial for the eigenvalue

$$\begin{align} \displaystyle (\lambda -1)\left(4\lambda ^3+5 \lambda ^2- (\mu +\sigma +3)\lambda-\mu \sigma +2\right) =0. \label{char2} \nonumber \end{align} \tag{ 41 }$$

The corresponding left boundary vectors are

$$\begin{align} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle &\bra{l_0}={\left(0,\frac{2 (\lambda_{\rm{L}} (\mu -1)-\mu) (\mu +\nu)}{\lambda_{\rm{L}} (2 \lambda_{\rm{L}}+\mu)}, \frac{\mu-\lambda_{\rm{L}}(\mu-1)}{\lambda_{\rm{L}}^2 (2 \lambda_{\rm{L}}-1)}, \frac{2 \lambda_{\rm{L}} (\nu +\mu (\mu +\nu)+1)-(\mu -1) \mu}{2 \lambda_{\rm{L}}+\mu } \right)}, \nonumber \\ &\bra{l_1}= {\left(0,\frac{2(\lambda_{\rm{L}} (\mu -1)-\mu) (\mu +\nu +1)}{\lambda_{\rm{L}}(2 \lambda_{\rm{L}}+\mu)}, -\frac{\mu-\lambda_{\rm{L}}(\mu-1)}{\lambda_{\rm{L}}^2 (2 \lambda_{\rm{L}}-1)}, \frac{(\mu +1) [(2 \lambda_{\rm{L}}-1) \mu +2 \lambda_{\rm{L}} \nu ]}{2 \lambda_{\rm{L}}+\mu } \right)}, \nonumber \\ &\bra{l^{\prime}_{0,0}}= \left(\mu -\frac{\mu }{\lambda_{\rm{L}}}-1, \mu ^2-2 (\lambda_{\rm L} (\mu-1) +1) (\mu +\nu)-1, 0, 0 \right), \nonumber \\ &\bra{l^{\prime}_{0,1}}= \left(0, 0, \lambda_{\rm{L}} (2 \lambda_{\rm{L}}-1) (\mu -1), \mu-\frac{\mu}{\lambda_{\rm{L}}} -1\right), \nonumber \\ &\bra{l^{\prime}_{1,0}}=\left(2 \lambda_{\rm{L}} (\mu -1)+\frac{\mu }{\lambda_{\rm{L}}}-3 \mu +1, 2 \lambda_{\rm{L}} (\mu -1) (\mu +\nu +1)-\mu ^2+2 \nu +1, 0, 0 \right), \nonumber \\ &\bra{l^{\prime}_{1,1}} =\left(0, 0, -\lambda_{\rm{L}} (\mu +1), \mu-\frac{\mu }{\lambda_{\rm{L}}} -1 \right), \nonumber \end{align} \tag{ 42 }$$

and the right boundary vectors are

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle &\ket{r_{0,0}}= \left(0, -\lambda_{\rm{R}} (\sigma +1), -2 \rho (\sigma \lambda_{\rm{R}}-\lambda_{\rm{R}}+1)-\sigma (2 \sigma \lambda_{\rm{R}}-2 \lambda_{\rm{R}}-\sigma +2)-1, 0 \right), \nonumber \\ &\ket{r_{0,1}}=\left(2 \lambda_{\rm{R}}^2 (\sigma +1), \lambda_{\rm{R}} (\sigma +1), -\sigma ^2+2\rho +2 \lambda_{\rm{R}} \left(\sigma ^2+\rho \sigma -\rho -1\right)+1, 0 \right), \nonumber \\ &\ket{r_{1,0}}=\left(0, \lambda_{\rm{R}} (2 \lambda_{\rm{R}}-1) (\sigma -1), 0, -\lambda_{\rm{R}}^2 (2 \lambda_{\rm{R}}-1)\right), \nonumber \\ &\ket{r_{1,1}} =\left(0,-\lambda_{\rm{R}} (\sigma +1), 0, -\lambda_{\rm{R}}^2 (2 \lambda_{\rm{R}}-1) \right), \nonumber \\ &\ket{r^{\prime}_0}=\left(-\lambda_{\rm{R}}^2,0, -\frac{\lambda_{\rm{R}} \left(2 \lambda_{\rm{R}} \left(\rho \sigma +\rho +\sigma ^2+1\right)-\sigma ^2+\sigma \right)}{2 \lambda_{\rm{R}}+\sigma }, \frac{2 \lambda_{\rm{R}}^3 (2\lambda_{\rm{R}}-1) (\rho +\sigma)}{2 \lambda_{\rm{R}}+\sigma }\right), \nonumber \\ &\ket{r^{\prime}_1}=\left(0, -\lambda_{\rm{R}}^2, \frac{\lambda_{\rm{R}} (\sigma +1) (2\lambda_{\rm{R}} (\rho +\sigma)-\sigma)}{2 \lambda_{\rm{R}}+\sigma }, -\frac{2\lambda_{\rm{R}}^3 (2 \lambda_{\rm{R}}-1) (\rho +\sigma +1)}{2\lambda_{\rm{R}}+\sigma } \right). \nonumber \end{align} \tag{ 43 }$$

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

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

$$\begin{align} \newcommand{\x}{{\rm x}} \displaystyle &F_{+}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & z \nonumber \\ 0 & 0 & \frac{\xi \omega -1}{\omega z^2} & 0 \nonumber \\ 0 & 0 & 0 & \xi z^2 \nonumber \\ \end{array} \right), \qquad \qquad \quad \; F_{-}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & \frac{1}{\xi ^2 z^3} \nonumber \\ 0 & 0 & \frac{\xi \omega -1}{\xi z^2} & 0 \nonumber \\ 0 & 0 & 0 & \omega +\frac{1}{\xi}\left(\frac{1}{z^2}-1\right) \nonumber \\ \end{array} \right), \nonumber \\ &F^{\prime}_{+}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & \frac{z^3}{\omega ^2} \nonumber \\ 0 & 0 & \frac{z^2 (\xi \omega -1)}{\omega } & 0 \nonumber \\ 0 & 0 & 0 & \xi +\frac{z^2-1}{\omega } \nonumber \\ \end{array} \right),\qquad F^{\prime}_{-}= \left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & \frac{1}{z} \nonumber \\ 0 & 0 & \frac{z^2 (\xi \omega -1)}{\xi } & 0 \nonumber \\ 0 & 0 & 0 & \frac{\omega }{z^2} \nonumber \\ \end{array} \right), \nonumber \\ &G_{+}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \xi ^2 & 0 & 0 & 0 \nonumber \\ -\frac{2 \xi }{z} & -1 & 0 & 0 \nonumber \\ \end{array} \right), \qquad \qquad \quad \;\, G_{-}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \xi \omega & 0 & 0 & 0 \nonumber \\ -\frac{\xi \omega z^2+1}{\xi z^3} & -\frac{\omega }{\xi z^2} & 0 & 0 \nonumber \\ \end{array} \right), \nonumber \\ &G^{\prime}_{+}= \left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \xi \omega & 0 & 0 & 0 \nonumber \\ -z^3 \left(\frac{\xi }{z^2}+\frac{1}{\omega}\right) & -\frac{\xi z^2}{\omega } & 0 & 0 \nonumber \\ \end{array} \right), \quad G^{\prime}_{-}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 0 \nonumber \\ \omega ^2 & 0 & 0 & 0 \nonumber \\ -2 \omega z & -1 & 0 & 0 \nonumber \\ \end{array} \right),\nonumber \\ &K_{+}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & \xi z^2 & 0 & 0 \nonumber \\ 0 & 0 & \frac{z^2 (\xi \omega -1)}{\omega } & 0 \nonumber \\ 0 & z & 0 & 0 \nonumber \\ \end{array} \right), \qquad \quad \;\;\; K_{-}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & \omega +\frac{1}{\xi}\left(\frac{1}{z^2}-1\right) & 0 & 0 \nonumber \\ 0 & 0 & \frac{z^2 (\xi \omega -1)}{\xi } & 0 \nonumber \\ 0 & \frac{1}{\xi ^2 z^3} & 0 & 0 \nonumber \\ \end{array} \right), \nonumber \\ &L_{+}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & \xi & 0 & -z \nonumber \\ \frac{\xi \omega }{z} & 0 & 0 & 0 \nonumber \\ 0 & 0 & -\frac{\xi \omega +z^2}{\omega ^2 z} & -\frac{z^2}{\omega } \nonumber \\ \end{array} \right),\qquad \quad \;\, L_{-}=\left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 0 \nonumber \\ 0 & \frac{1}{\xi z^2} & 0 & -\frac{\omega }{\xi z} \nonumber \\ 0 & 0 & \frac{\omega }{z^2} & 0 \nonumber \\ -2 \omega & 0 & 0 & -\omega \nonumber \\ \end{array} \right). \nonumber \end{align} \tag{ 44 }$$

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

$$\begin{align} \newcommand{\x}{{\rm x}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle &\hat{W}_s= e_{11}\otimes W_s(\xi z, \omega/z)+ e_{22}\otimes W_s(\xi/ z, \omega z), \nonumber \\ & \hat{W}^{\prime}_s=e_{11}\otimes W^{\prime}_s(\xi z, \omega/z)+ e_{22}\otimes W^{\prime}_s(\xi/ z, \omega z), \label{diagpart} \nonumber \end{align} \tag{ 45 }$$

where $\renewcommand{\bra}[1]{{\langle #1 \vert}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} e_{i j} = \ket{i}\bra{\,j}$, $i, j\in \{1, 2\}$, denotes a basis of $2\times 2$ matrices. Further we define the following block triangular $8\times 8$ matrices in extended auxiliary space

$$\begin{align} \newcommand{\one}{\mathbb {1}} \newcommand{\x}{{\rm x}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle &\hat{F}^{(k)}=\one_8+ e_{1 2} \otimes \frac{c_+ z^k F_+ +c_- z^{-k} F_-}{\xi \omega -1}, \quad \hat{F}^{\prime (k)}=\one_8+ e_{1 2} \otimes \frac{c_+ z^k F^{\prime}_+ +c_- z^{-k} F^{\prime}_-}{\xi \omega -1}, \nonumber \\ &\hat{G}^{(k)}=\one_8+ e_{1 2} \otimes \frac{c_+ z^k G_+ +c_- z^{-k} G_-}{\xi \omega -1}, \quad \hat{G}^{\prime (k)}=\one_8+ e_{1 2} \otimes \frac{c_+ z^k G^{\prime}_+ +c_- z^{-k} G^{\prime}_-}{\xi \omega -1}, \nonumber \\ &\hat{K}^{(k)}=\one_8+ e_{1 2} \otimes \frac{c_+ z^k K_+ +c_- z^{-k} K_-}{\xi \omega -1}, \nonumber \\ &\hat{L}^{(k)}=(z e_{11}+z^{-1} e_{22})\otimes \one_4+ e_{1 2} \otimes \frac{c_+ z^k L_+ +c_- z^{-k} L_-}{\xi \omega -1}. \label{FLhat} \nonumber \end{align} \tag{ 46 }$$

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

$$\begin{align} \displaystyle &P_{123} \hat{K}^{(k-1)}\hat{\bf {W}}_1 \hat{S} \hat{G}^{(k)}\hat{\bf {W}}_2 \hat{F}^{(k+1)}\hat{\bf {W}}^{\prime}_3=\hat{F}^{\prime (k-1)}\hat{\bf {W}}_1 \hat{G}^{\prime (k)} \hat{\bf {W}}^{\prime}_2 \hat{K}^{(k+1)}\hat{\bf {W}}_3 \hat{S} \nonumber \\ &P_{123} \hat{G}^{\prime (k-1)}\hat{\bf {W}}^{\prime}_1 \hat{F}^{\prime (k)} \hat{\bf {W}}_2 \hat{L}^{(k+1)}\hat{\bf {W}}^{\prime}_3 \hat{S}=\hat{L}^{(k-1)}\hat{\bf {W}}^{\prime}_1 \hat{S} \hat{F}^{(k)}\hat{\bf {W}}^{\prime}_2 \hat{G}^{(k+1)}\hat{\bf {W}}_3, \label{bulk2} \nonumber \end{align} \tag{ 47 }$$

where $\newcommand{\ot}[1]{{\widetilde{#1}}} \newcommand{\one}{\mathbb {1}} \hat{S}=\one_2 \otimes S$, which can be straightforwardly checked to hold for any $\newcommand{\CC}{\mathbb{C}} \newcommand{\x}{{\rm x}} \xi, \omega, z, c_+, c_-\in\CC$ and $\newcommand{\ZZ}{\mathbb{Z}} k\in\ZZ$. Note that by setting $z=1, c_+=c_-=0$, we recover the original bulk algebra (18) and (20).

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

$$\begin{align} \displaystyle \mathcal{R} \hat{F}^{(k)}=\hat{F}^{\prime (k)} \qquad \mathcal{R} \hat{G}^{(k)}=\hat{G}^{\prime (k)} \qquad \mathcal{R} \hat{{\bf W}}=\hat{{\bf W}}^{\prime}. \nonumber \end{align} \tag{ 48 }$$

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

Lemma. Let us assume that 8-dimensional boundary vectors $\renewcommand{\bra}[1]{{\langle #1 \vert}} \bra{\hat{l}_{s}}$, $\renewcommand{\bra}[1]{{\langle #1 \vert}} \bra{\hat{l}^{\prime}_{s, s^{\prime}}}$, $\renewcommand{\ket}[1]{{\vert #1 \rangle}} \ket{\hat{r}_{s, s^{\prime}}}$, $\renewcommand{\ket}[1]{{\vert #1 \rangle}} \ket{\hat{r}^{\prime}_s}$ exist, together with parameters $\newcommand{\x}{{\rm x}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \xi, \omega, z, c_+, c_-, \Lambda_{\rm L}, \Lambda_{\rm R}$, such that the following boundary equations are satisfied

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hmp}[2]{{\hat{#1}^{^{\prime}(#2)}}} \newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle P^{\rm{L}}_{12} \bra{\hat{\bf {l}}^{\prime}_{12}}=\Lambda_{\rm{L}} \bra{\hat{\bf {l}}_1} \hmp{G}0 \hmup{W}{2}, \label{bound1} \nonumber \end{align} \tag{ 49 }$$

$$\begin{align} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle P_{123} \bra{\hat{\bf {l}}_{1}} \hm{L}{0}\hmup{W}{2}\hat{S} \hm{F}{1}\hmup{W}{3}=\bra{\hat{\bf {l}}^{\prime}_{12}} \hm{K}{1}\hmu{W}{3}\hat{S},\label{bound2} \nonumber \end{align} \tag{ 50 }$$

$$\begin{align} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle P^{R}_{12} \ket{\hat{\bf {r}}_{12}}=\hm{F}{n-3}\hmup{W}{1}\hat{S}\ket{\hat{\bf {r}}^{\prime}_2}, \label{bound3} \nonumber \end{align} \tag{ 51 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle P_{123} \hm{G}{n-4}\hmup{W}{1} \hm{K}{n-3}\hmu{W}{2}\hat{S} \ket{\hat{\bf {r}}^{\prime}_{3}} =\Lambda_{\rm{R}} \hm{L}{n-4}\hmup{W}{1}\hat{S}\ket{\hat{\bf {r}}_{23}}. \label{bound4} \nonumber \end{align} \tag{ 52 }$$

Then, the following inhomogeneous (site-dependent) MPA

$$\begin{align} \newcommand{\WW}{{\hat{{\bf W}}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p}=\bra{\hat{\bf {l}}_1}\hat{L}^{(0)}\WW^{\prime}_2\hat{S}\hat{F}^{(1)}{\WW}^{\prime}_3\hat{G}^{(2)}{\WW}_4 \cdots \hat{F}^{(n-5)}{\WW}^{\prime}_{n-3}\hat{G}^{(n-4)}{\WW}_{n-2}\ket{\hat{\bf {r}}_{n-1,n}}, \label{p} \nonumber \end{align} \tag{ 53 }$$

$$\begin{align} \newcommand{\WW}{{\hat{{\bf W}}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p^{\prime}}=\bra{\hat{\bf {l}}^{\prime}_{12}}\hat{F}^{^{\prime}(1)}{\WW}_3\hat{G}^{^{\prime}(2)}{\WW}^{\prime}_4\hat{F}^{^{\prime}(3)}{\WW}_5 \cdots \hat{G}^{^{\prime}(n-4)}{\WW}^{\prime}_{n-2}\hat{K}^{(n-3)}{\WW}_{n-1}\hat{S}\ket{\hat{\bf {r}}^{\prime}_{n}}, \label{pp} \nonumber \end{align} \tag{ 54 }$$

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

Proof. Explaining how cancellation mechanism works is fully analogous to the simpler case of spatially homogeneous NESS-orbital (21)–(26). However, due to commutativity (9) we can explain it here for the reverse order: when U_e acts on ${\bf p}$, $P_{123}$ first acts on the left boundary vector and via (50) creates $\newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \hm{K}{1}\hmu{W}{3}\hat{S}$. The subsequent P's in U_e transfer $\newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \hm{K}{1}\hmu{W}{3}\hat{S}$ to the right via the bulk algebra (47). Before the final $P_{n-3, n-2, n-1}$ acts, $P^{\rm{R}}_{n-1, n}$ acts (as it commutes with $P_{n-3, n-2, n-1}$) and creates $\newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \hm{F}{n-3}\hmup{W}{1}\hat{S}$ necessary for the final $P_{n-3, n-2, n-1}$ to transfer $\newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \hm{K}{n-5}\hmu{W}{n-3}\hat{S}$ to the right end as $\newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \hm{K}{n-3}\hmu{W}{n-1}\hat{S}$, thus finally creating ${\bf p^{\prime}}$ (54). The odd part of the propagator U_o acts analogously in reverse. □

We illustrate this mechanism by writing out an explicit example for $n=8$ and for the even part of the propagator U_e,

$$\begin{align} \newcommand{\WW}{{\hat{{\bf W}}}} \newcommand{\hmup}[2]{{\hat{{\bf #1}}^{\prime}_{#2}}} \newcommand{\hmu}[2]{{\hat{{\bf #1}}_{#2}}} \newcommand{\hm}[2]{{\hat{#1}^{(#2)}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle U_{\rm e} {\bf p} &= P_{123}P_{345}P_{567} P^{\rm{R}}_{78} \bra{\hat{\bf {l}}_1}\hat{L}^{(0)}\WW^{\prime}_2\hat{S}\hat{F}^{(1)}{\WW}^{\prime}_3\hat{G}^{(2)}{\WW}_4 \hat{F}^{(3)}{\WW}^{\prime}_{5}\hat{G}^{(4)}{\WW}_{6}\ket{\hat{\bf {r}}_{78}}\nonumber \\ &=P_{345}P_{567} P^{\rm{R}}_{78}\bra{\hat{\bf {l}}^{\prime}_{12}} \hm{K}{1}\hmu{W}{3}\hat{S}\hat{G}^{(2)}{\WW}_4 \hat{F}^{(3)}{\WW}^{\prime}_{5}\hat{G}^{(4)}{\WW}_{6}\ket{\hat{\bf {r}}_{78}}\nonumber \\ &=P_{567} P^{\rm{R}}_{78}\bra{\hat{\bf {l}}^{\prime}_{12}} \hat{F}^{\prime (1)}\hat{\bf {W}}_1 \hat{G}^{\prime (2)} \hat{\bf {W}}^{\prime}_2 \hat{K}^{(3)}\hat{\bf {W}}_3 \hat{S}\hat{G}^{(4)}{\WW}_{6}\ket{\hat{\bf {r}}_{78}}\nonumber \\ &=P_{567} \bra{\hat{\bf {l}}^{\prime}_{12}} \hat{F}^{\prime (1)}\hat{\bf {W}}_1 \hat{G}^{\prime (2)} \hat{\bf {W}}^{\prime}_2 \hat{K}^{(3)}\hat{\bf {W}}_3 \hat{S}\hat{G}^{(4)}{\WW}_{6}\hm{F}{5}\hmup{W}{7}\hat{S}\ket{\hat{\bf {r}}^{\prime}_8}\nonumber \\ &=\bra{\hat{\bf {l}}^{\prime}_{12}} \hat{F}^{\prime (1)}\hat{\bf {W}}_1 \hat{G}^{\prime (2)} \hat{\bf {W}}^{\prime}_2 \hat{F}^{\prime (3)}\hat{\bf {W}}_1 \hat{G}^{\prime (4)} \hat{\bf {W}}^{\prime}_2 \hat{K}^{(5)}\hat{\bf {W}}_3 \hat{S}^2\ket{\hat{\bf {r}}^{\prime}_8} \nonumber \\ & = {\bf p^{\prime}}. \nonumber \end{align} \tag{ 55 }$$

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

To solve the boundary equations we follow a similar procedure as outlined in section 2 for the NESS-orbital. Let us decompose the extended auxiliary space as a direct sum of two 4-dimensional spaces ${{\mathcal H}}_1\oplus {{\mathcal H}}_2$, where an element of ${{\mathcal H}}_i$ is written as $\newcommand{\ot}[1]{{\widetilde{#1}}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \ket{i}\otimes \ket{\psi}$. Due to the upper triangular structure (46) the left boundary equations projected to the subspace ${{\mathcal H}}_1$ reduce to those for the NESS-orbital with a scaling factor z coming from $e_{11}$ component of $\hat{L}^{(k)}$, but with rescaled spectral parameters ($\newcommand{\x}{{\rm x}} \xi\to \xi z$, $\omega \to \omega/z$, because of (45)). Since the NESS solution we found in section 2 is unique, the left boundary vector in this subspace must be the NESS-orbital boundary vector with scaled $\newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \lambda_{\rm L} = \Lambda_{\rm L}/z$. Comparing (49), (50) with (23), (24), (28) and (29) immediately fixes the values of the spectral parameters ξ and ω (writing first for the conditional driving (13)):

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{z (\alpha +\beta -1)-\beta \Lambda_{\rm{L}}}{(\beta -1) \Lambda_{L}^{2}}, \label{xi1} \nonumber \end{align} \tag{ 56 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \omega =\frac{\Lambda_{\rm{L}} (\alpha z-\Lambda_{\rm{L}})}{(\beta -1) z}. \label{paraleft1} \nonumber \end{align} \tag{ 57 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{c_-}{c_+}=\frac{\Lambda_{L}^{4}}{z^4 (\alpha +\beta -1)}, \label{paraleft2} \nonumber \end{align} \tag{ 58 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle &\bra{\hat{l}_s}=\bra{1} \otimes \Big\langle l_s\Big(\lambda_{\rm{L}}=\frac{\Lambda_{\rm{L}}}{z}\Big)\Big| + c_+ \bra{2} \otimes \bra{\tilde{l}_s},\nonumber \\ &\bra{\hat{l}^{\prime}_{s,s^{\prime}}}=z \bra{1} \otimes \Big\langle l^{\prime}_{s,s^{\prime}}\Big(\lambda_{\rm{L}}=\frac{\Lambda_{\rm{L}}}{z}\Big)\Big|+ c_+ \bra{2} \otimes \bra{\tilde{l}^{\prime}_{s,s^{\prime}}}, \label{leftvecs} \nonumber \end{align} \tag{ 59 }$$

where the explicit form of the 'offdiagonal' vectors $\renewcommand{\bra}[1]{{\langle #1 \vert}} \bra{\tilde{l}_s}$ and $\renewcommand{\bra}[1]{{\langle #1 \vert}} \bra{\tilde{l}^{\prime}_{s, s^{\prime}}}$ are given in appendix A.1. Following a fully analogous procedure for the right boundary equations (51) and (52), we arrive to the following results for the spectral parameters and the amplitude ratio

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\Lambda_{\rm{R}} (\gamma z-\Lambda_{\rm{R}})}{(\delta -1) z}, \label{xi2} \nonumber \end{align} \tag{ 60 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \omega =\frac{z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}}}{(\delta -1) \Lambda_{\rm{R}}^2},\label{om2} \nonumber \end{align} \tag{ 61 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{c_+}{c_-}=\frac{\Lambda_{\rm{R}}^{4} z^{4 m+2}}{\gamma +\delta -1}, \label{pararight} \nonumber \end{align} \tag{ 62 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \newcommand{\ot}[1]{{\widetilde{#1}}} \displaystyle &\ket{\hat{r}_{s,s^{\prime}}}= c_- z^{n-3}\ket{1} \otimes \ket{\tilde{r}_{s,s^{\prime}}} + \ket{2} \otimes \Big| r_{s,s^{\prime}}\Big(\lambda_{\rm{R}}=\frac{\Lambda_{\rm{R}}}{z} \Big)\Big\rangle, \nonumber \\ &\ket{\hat{r}^{\prime}_s}=c_- z^{n-3} \ket{1} \otimes \ket{\tilde{r}^{\prime}_s} + \ket{2} \otimes \Big| r^{\prime}_s\Big(\lambda_{\rm{R}}=\frac{\Lambda_{\rm{R}}}{z}\Big)\Big\rangle, \label{rightvecs} \nonumber \end{align} \tag{ 63 }$$

whose components in ${{\mathcal H}}_2$ are expressed in terms of right boundary vectors for the NESS-orbital and the complementary (offdiagonal) components $\renewcommand{\ket}[1]{{\vert #1 \rangle}} \ket{\tilde{r}_{s, s^{\prime}}}$ and $\renewcommand{\ket}[1]{{\vert #1 \rangle}} \ket{\tilde{r}^{\prime}_{s}}$ are given in appendix A.1. Pairwise identifying equations (56), (60), (57), (61), (58) and (62) we obtain a closed set of Bethe-like equations for $\newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \Lambda_{\rm L}, \Lambda_{\rm R}, z$:

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{z (\alpha +\beta -1)-\beta \Lambda_{\rm{L}}}{(\beta -1) \Lambda_{L}^{2}}=\frac{\Lambda_{\rm{R}} (\gamma z-\Lambda_{\rm{R}})}{(\delta -1) z}, \nonumber \end{align} \tag{ 64 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{L}} (\alpha z-\Lambda_{\rm{L}})}{(\beta -1) z}=\frac{z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}}}{(\delta -1) \Lambda_{\rm{R}}^2}, \nonumber \end{align} \tag{ 65 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle (\alpha +\beta -1)(\gamma +\delta -1)=\Lambda_{\rm{L}}^4 \Lambda_{\rm{R}}^4 z^{4m-2}. \nonumber \end{align} \tag{ 66 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle &\Lambda ^4-\alpha \gamma \Lambda ^3 z^2+ (\alpha +\beta +\gamma +\delta-\alpha\delta-\beta \gamma-2)\Lambda ^2 z^2 \nonumber \\ &\quad -\beta \delta \Lambda z^2+ (\alpha +\beta -1) (\gamma +\delta -1)z^4=0,\label{bethe1} \nonumber \end{align} \tag{ 67 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle (\alpha +\beta -1)(\gamma +\delta -1)-\Lambda^4 z^{4m-2}=0. \label{bethe2} \nonumber \end{align} \tag{ 68 }$$

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

For the case of Bernoulli driving (14) we find the following solutions for the left boundary equations (49) and (50)

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\Lambda_{\rm{L}}+\mu \left(z-\Lambda_{\rm{L}}\right)}{(\mu +1) \Lambda_{\rm{L}}^2}, \nonumber \end{align} \tag{ 69 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \omega=\frac{\Lambda_{\rm{L}} \left(2 \Lambda_{\rm{L}}-z\right)}{z (\mu +1)}, \nonumber \end{align} \tag{ 70 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{c_-}{c_+}=\frac{2 \Lambda_{\rm{L}}^4}{z^4 \mu}, \nonumber \end{align} \tag{ 71 }$$

and for the right ones (51) and (52)

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \newcommand{\x}{{\rm x}} \displaystyle \xi=\frac{\Lambda_{\rm{R}} \left(2 \Lambda _{\rm{R}}-z \right)}{z (\sigma+1)}, \nonumber \end{align} \tag{ 72 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \omega=\frac{\Lambda_{\rm{R}}+\sigma \left(z-\Lambda_{\rm{R}}\right)}{(\sigma +1) \Lambda_{\rm{R}}^2}, \nonumber \end{align} \tag{ 73 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{c_+}{c_- }=\frac{2 z^{4m+2} \Lambda_{\rm{R}}^4}{\sigma }. \nonumber \end{align} \tag{ 74 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{L}}+\mu \left(z-\Lambda_{\rm{L}}\right)}{(\mu +1) \Lambda_{\rm{L}}^2}=\frac{\Lambda_{\rm{R}} \left(2 \Lambda_{\rm{R}} -z\right)}{z (\sigma +1)} \nonumber \end{align} \tag{ 75 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{R}}+\sigma \left(z-\Lambda_{\rm{R}}\right)}{(\sigma +1) \Lambda_{\rm{R}}^2}=\frac{\Lambda_{\rm{L}} \left(2 \Lambda_{\rm{L}}-z\right)}{z (\mu +1)} \nonumber \end{align} \tag{ 76 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{1}{4} \mu \sigma= \Lambda_{\rm{L}}^4 \Lambda_{\rm{R}}^4 z^{4 m-2}, \nonumber \end{align} \tag{ 77 }$$

or equivalently

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle 4 \Lambda ^4+\Lambda ^3 z^2+(\mu+\sigma -2)\Lambda ^2 z^2- (\mu -1) (\sigma -1)\Lambda z^2+\mu \sigma z^4 =0 \label{bethe1a}, \nonumber \end{align} \tag{ 78 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{1}{4} \mu \sigma- \Lambda^4 z^{4 m-2}=0.\label{bethe2a} \nonumber \end{align} \tag{ 79 }$$

The boundary vectors in the Bernoulli driving case are of the same form as in the conditional driving case, namely (59) and (63). The explicit expressions for their components are given in appendix A.2.

3.1. Thermodynamics

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

$$\begin{align} \newcommand{\ii}{{\rm i}} \displaystyle z={\rm e}^{\ii \kappa}, \label{uniz} \nonumber \end{align} \tag{ 80 }$$

where κ is the quasi-momentum which may be restricted to interval $[0, \pi)$ as (67) and (68) (or (78) and (79)) are symmetric under the transformation $z \to -z$. Thus, in the leading order of $1/n$ expansion the spectrum in the first orbital converges to an algebraic curve $\newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \Lambda(\kappa)$ given by (80) and (67) (or (78)) (see figure 3).

The case $\newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \Lambda(\kappa=0)=1$ corresponds to the eigenvalues which in the thermodynamic limit converge to 1 and their scaling with $1/n$ will determine the asymptotic relaxation rate of the system. Writing the expansion of $\newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\L}{{\rm{L}}} \Lambda(0)$ as

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \Lambda(0)=1-\Lambda_1(0)/n+\mathcal{O}(n^{-2}) \nonumber \end{align} \tag{ 81 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \Lambda_1(0)=\frac{[1-(\alpha +\beta -1) (\gamma +\delta -1)] \log [(\alpha +\beta -1) (\gamma+\delta -1)]}{2 (\alpha +\beta -1) (\gamma +\delta -1)+\alpha \gamma -\beta \delta-2}, \label{ldmth} \nonumber \end{align} \tag{ 82 }$$

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \Lambda_1(0)=\frac{[4-\mu \sigma] \log \left[\frac{1}{4} \mu \sigma \right]}{(\mu +1) (\sigma +1)-9}, \label{ldmtha} \nonumber \end{align} \tag{ 83 }$$

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

3.2. Quadratic form of the bulk algebra

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

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

$$\begin{align} \displaystyle &\hat{{\bf Z}}_{12}^{(k)}=\hat{F}^{(2k-1)} \hat{{\bf W}}^{\prime}_1 \hat{G}^{(2k)} \hat{{\bf W}}_2, \nonumber \\ &\hat{{\bf Z}}_{12}^{^{\prime}(k)}=\hat{F}^{^{\prime}(2k-1)} \hat{{\bf W}}_1 \hat{G}^{^{\prime}(2k)} \hat{{\bf W}}^{\prime}_2, \nonumber \\ &\hat{{\bf Y}}_{12}^{(k)}=\hat{K}^{(2k-1)} \hat{{\bf W}}_1 \hat{S} \hat{G}^{(2k)} \hat{{\bf W}}_2, \nonumber \\ &\hat{{\bf Y}}_{12}^{^{\prime}(k)}=\hat{F}^{^{\prime}(2k-1)} \hat{{\bf W}}_1 \hat{L}^{(2k)} \hat{{\bf W}}^{\prime}_2 \hat{S}. \nonumber \end{align} \tag{ 84 }$$

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

$$\begin{align} \displaystyle &P_{123}\hat{{\bf Y}}_{12}^{(k)} \hat{{\bf Z}}_{34}^{(k+1)}=\hat{{\bf Z}}_{12}^{^{\prime}(k)}\hat{{\bf Y}}_{34}^{(k+1)}, \nonumber \\ &P_{234} \hat{{\bf Z}}_{12}^{^{\prime}(k)}\hat{{\bf Y}}_{34}^{^{\prime}(k+1)}=\hat{{\bf Y}}_{12}^{^{\prime}(k)} \hat{{\bf Z}}_{34}^{(k+1)}. \nonumber \end{align} \tag{ 85 }$$

This quadratic formulation has perhaps some appeal as it is reminiscent of ZF algebra⁸. The eigenvector MPA (53) and (54) now transforms to:

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p}=\bra{\hat{\bf {l}}_{12}}\hat{{\bf Z}}_{34}^{(1)}\hat{{\bf Z}}_{56}^{(2)} \ldots \hat{{\bf Z}}_{n-3,n-2}^{(m)} \ket{\hat{\bf {r}}_{n-1,n}}, \nonumber \end{align} \tag{ 86 }$$

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle {\bf p^{\prime}}=\bra{\hat{\bf {l}}^{\prime}_{12}}\hat{{\bf Z}}_{34}^{^{\prime}(1)}\hat{{\bf Z}}_{56}^{^{\prime}(2)} \ldots \hat{{\bf Z}}_{n-3,n-2}^{^{\prime}(m)}\ket{\hat{\bf {r}}^{\prime}_{n-1,n}}, \nonumber \end{align} \tag{ 87 }$$

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

$$\begin{align} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\hat{\bf {l}}_{12}}= \bra{\hat{\bf {l}}_{1}} \hat{L}^{(0)} \hat{{\bf W}}^{\prime}_2 \hat{S}, \nonumber \end{align} \tag{ 88 }$$

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\hat{\bf {r}}^{\prime}_{n-1,n}}=\hat{K}^{(n-3)} \hat{{\bf W}}_{n-1} \hat{S} \ket{\hat{\bf {r}}^{\prime}_{n}}. \nonumber \end{align} \tag{ 89 }$$

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

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

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

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{z (\alpha +\beta -1)-\beta \Lambda_{\rm{L}}}{(\beta -1) \Lambda_{L}^{2}}=\frac{\Lambda_{\rm{R}} (\gamma z-\Lambda_{\rm{R}})}{(\delta -1) z}, \nonumber \end{align} \tag{ 90 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{L}} (\alpha z-\Lambda_{\rm{L}})}{(\beta -1) z}=\frac{z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}}}{(\delta -1) \Lambda_{\rm{R}}^2}, \nonumber \end{align} \tag{ 91 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle (\alpha +\beta -1)^{\,p}(\gamma +\delta -1)^{\,p}=(\Lambda_{\rm{L}} \Lambda_{\rm{R}})^{4 p} z^{4(m-p)+2}, \nonumber \end{align} \tag{ 92 }$$

or for the Bernouli driving (14):

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{L}}+\mu \left(z-\Lambda_{\rm{L}}\right)}{(\alpha -\beta +1) \Lambda_{\rm{L}}^2}=\frac{\Lambda_{\rm{R}} \left(2 \Lambda_{\rm{R}} -z\right)}{z (\sigma +1)}, \nonumber \end{align} \tag{ 93 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \frac{\Lambda_{\rm{R}}+\sigma \left(z-\Lambda_{\rm{R}}\right)}{(\sigma +1) \Lambda_{\rm{R}}^2}=\frac{\Lambda_{\rm{L}} \left(2 \Lambda_{\rm{L}}-z\right)}{z (\mu +1)}, \nonumber \end{align} \tag{ 94 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle 4^{-p} \mu^{\,p} \sigma^{\,p}=(\Lambda_{\rm{L}} \Lambda_{\rm{R}})^{4 p} z^{4(m-p)+2}. \nonumber \end{align} \tag{ 95 }$$

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

Furthermore, we provide a conjecture for the average degeneracy of the eigenvalues of the pth orbital which seems to scale as

$$\begin{align} \displaystyle g(m,p)=\frac{1}{p} \left(\begin{array}{@{}c@{}} 2 m \nonumber \\ p-1 \end{array}\right). \nonumber \end{align} \tag{ 96 }$$

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

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

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

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

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

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

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

In this appendix we list the so-called off-diagonal components of the boundary vectors which solve the boundary equations (49)–(52). We note that the values of all other parameters (which enter into the equations in a nonlinear way) and other ('diagonal') components of the boundary vectors can be fixed by comparing to a $4$-dimensional auxiliary problem for the NESS-orbital. The remaining components reported below are then fixed by solving the remaining system of linear equations. They are quite lengthy and their algebraic form could probably still be considerably optimized, but they are given here just for completeness.

A.1. Boundary vectors for conditional driving

The offdiagonal left boundary vector components for the first orbital (59) and for conditional driving (13) are given as

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle &\bra{\tilde{l}_{0}}= C_1 \bigg(0,\frac{(\beta -1)^2 \beta \Lambda_{\rm{L}} z \left(\Lambda_{\rm{L}}-z^3 (\alpha+\beta -1)\right)}{\left(z^2-1\right) (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}, \nonumber \\ &c_1, -\frac{(\beta -1)^2 \beta \Lambda_{\rm{L}}^2 z^2 \left(\beta \Lambda_{\rm{L}}^2+\alpha ^2 z^2-\alpha \Lambda_{\rm{L}} z \left((\beta -1) z^2+2\right)\right)}{\left(z^2-1\right) (\alpha z-\Lambda_{\rm{L}}) (\alpha z-\beta \Lambda_{\rm{L}})} \bigg),\nonumber \end{align} \tag{ A.1 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle &\bra{\tilde{l}_{1}} = C_1\bigg(0, \frac{(\beta -1)^3 \Lambda_{\rm{L}} z \left(z^3 (\alpha +\beta -1)-\Lambda_{\rm{L}}\right)}{\left(z^2-1\right) (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}, \nonumber \\ &c_2, \frac{(\beta -1)^3 \Lambda_{\rm{L}}^2 z^2 \left(\beta \Lambda_{\rm{L}}^2+\alpha ^2 z^2-\alpha \Lambda_{\rm{L}} z \left((\beta -1) z^2+2\right)\right)}{\left(z^2-1\right) (\alpha z-\Lambda_{\rm{L}}) (\alpha z-\beta \Lambda_{\rm{L}})} \bigg), \nonumber \end{align} \tag{ A.2 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{0,0}} =\left(\frac{\alpha \Lambda_{\rm{L}} (\Lambda_{\rm{L}}+z (1-\alpha -\beta))}{\beta z \left(z^2-1\right) (\alpha +\beta -1)}, \frac{\Lambda_{\rm{L}}^2 (\Lambda_{\rm{L}}+z (1-\alpha -\beta))}{\left(z^2-1\right) (\alpha +\beta -1)}, 0, 0 \right), \nonumber \end{align} \tag{ A.3 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{0,1}}=\left(0, 0, c_3, \frac{\alpha \Lambda_{\rm{L}}(\Lambda_{\rm{L}}+z (1-\alpha -\beta))}{\beta z \left(z^2-1\right) (\alpha+\beta -1)} \right), \nonumber \end{align} \tag{ A.4 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{1,0}}= \left(\frac{\Lambda_{\rm{L}} (\alpha -\Lambda_{\rm{L}} z) (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}{\beta z \left(z^2-1\right) (\alpha +\beta -1)}, \frac{(\beta -1) \Lambda_{\rm{L}}^2 (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}{\beta \left(z^2-1\right) (\alpha +\beta -1)}, 0, 0 \right), \nonumber \end{align} \tag{ A.5 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{1,1}} =\left(0, 0, c_4, \frac{(\alpha -1)\Lambda_{\rm{L}} (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}{\beta z \left(z^2-1\right) (\alpha +\beta -1)}\right), \nonumber \end{align} \tag{ A.6 }$$

where

$$\begin{align*} \displaystyle C_1=\frac{1}{(\beta-1) (\alpha + \beta-1) \beta}, \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_1= \,\,& \beta -\frac{\alpha ^2 (\beta -1) z^4}{(\Lambda_{\rm{L}}-\alpha z)^2}+\frac{\alpha z \left(\alpha +\alpha ^2 z^2-\alpha z^2-\beta ^2 z^2\right)}{\beta (\alpha +\beta) (\beta \Lambda_{\rm{L}}-\alpha z)}-1\nonumber \\ &+\frac{\alpha }{\beta(1 - z^2)}-\frac{z(\alpha +\beta -1)^2 \left(\alpha +\beta +z^2 (\alpha +\beta +1)-1\right)}{(\alpha+\beta) (z (\alpha +\beta -1)-\Lambda_{\rm{L}})}\nonumber \\ &+\frac{\alpha z \left(-\alpha \beta +\alpha +\beta +\beta z^4 (\alpha +2 \beta-3)+z^2 \left(-\alpha -2 \beta ^2+\beta +2\right)-1\right)}{\beta \left(z^2-1\right)(\alpha z-\Lambda_{\rm{L}})}, \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_2=\,\,&\frac{\alpha ^2 (\beta -1)^2 z^4}{(\alpha z-\Lambda_{\rm{L}})^2}+\frac{z \left(\alpha^2-\beta +z^2 \left(\alpha ^3-(\alpha +1) \beta ^2+\beta \right)\right)}{(\alpha+\beta) (\alpha z-\beta \Lambda_{\rm{L}})}\nonumber \\ &+\frac{(\beta -1) z (\alpha +\beta-1)^2 \left(\alpha +\beta +z^2 (\alpha +\beta +1)-1\right)}{(\alpha +\beta) (z(\alpha +\beta -1)-\Lambda_{\rm{L}})}\nonumber \\ &+\frac{z \left[(\alpha -1) (\alpha (\beta -2)+1)-\alpha (\beta -1) z^4 (\alpha +2 \beta -3)\right]}{\left(z^2-1\right) (\alpha z-\Lambda_{\rm{L}})}\nonumber \\ &+\frac{z\left[z^2 (\alpha (\alpha +\beta (2 \beta-3)-1)-\beta +2)\right]}{\left(z^2-1\right) (\alpha z-\Lambda_{\rm{L}})}-\frac{\alpha +(\beta -2) \beta +(\beta -1)^2 z^2}{z^2-1}, \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_3=-\frac{\Lambda_{\rm{L}}^3 (\alpha z-\Lambda_{\rm{L}}) \left(\beta \Lambda_{\rm{L}}^2+z^4 (\alpha +\beta -1)^2-\Lambda_{\rm{L}} z (\alpha +\beta -1)\left(\beta +z^2\right)\right)}{z \left(z^2-1\right) (\alpha +\beta -1) (z (\alpha+\beta -1)-\beta \Lambda_{\rm{L}})^2}, \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_4=\,& \frac{\Lambda_{\rm{L}}^3(\alpha -1) (\beta -1) (z (\alpha +\beta-1)-\Lambda_{\rm{L}}) }{\beta \left(z^2-1\right) (\alpha +\beta -1) (\alpha z-\Lambda_{\rm{L}}) (z (\alpha +\beta -1)-\beta \Lambda_{\rm{L}})^2} \nonumber \\ &\times \left[\alpha z^4 (\alpha +\beta -1)-\beta \Lambda_{\rm{L}}^2 \left(z^2-2\right)-\Lambda_{\rm{L}} z (\alpha \beta +\alpha +\beta-1)\right]. \nonumber \end{align*}$$

The offdiagonal right boundary vectors components (63) are

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0,0}}=C_2\left(\frac{\gamma (1-\delta) \Lambda_{\rm{R}}}{\delta z \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, 0, c_5, 0 \right), \nonumber \end{align} \tag{ A.7 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0,1}}=C_2 \left(\frac{(\delta -1) \Lambda_{\rm{R}} (\gamma -\Lambda_{\rm{R}} z) }{\delta z \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, 0, c_6, 0\right), \nonumber \end{align} \tag{ A.8 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle &\ket{\tilde{r}_{1,0}}= C_2 \bigg(0, \frac{\Lambda_{\rm{R}} (\Lambda_{\rm{R}} z-\gamma)}{z \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, 0, \nonumber \\ &\frac{\gamma \Lambda_{\rm{R}}^3(\Lambda_{\rm{R}} z-\gamma) (z (\gamma +\delta -1)-\Lambda_{\rm{R}})}{\delta \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})} \bigg),\nonumber \end{align} \tag{ A.9 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle &\ket{\tilde{r}_{1,1}}= C_2\bigg(0, \frac{(\gamma -1) (\delta -1) \Lambda_{\rm{R}}}{\delta z \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, 0, \nonumber \\ &\frac{(\gamma -1) \Lambda_{\rm{R}}^3(\gamma -\Lambda_{\rm{R}} z) (z (\gamma +\delta -1)-\Lambda_{\rm{R}})}{\delta \left(z^2-1\right) (\gamma +\delta -1) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})} \bigg),\nonumber \end{align} \tag{ A.10 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0}}= C_3 \left(0, c_7, (\delta -1) \delta \Lambda_{\rm{R}}^2 \theta, (\delta -1) \delta \Lambda_{\rm{R}}^4 z^2 (\Lambda_{\rm{R}}-\gamma z)^2 (\delta \Lambda_{\rm{R}} z-\gamma) \right), \nonumber \end{align} \tag{ A.11 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{1}} =C_3 \left(0, \frac{1-\gamma}{\gamma} c_7,-(\delta -1)^2 \Lambda_{\rm{R}}^2 \theta, (\delta -1)^2 \Lambda_{\rm{R}}^4 z^2 (\Lambda_{\rm{R}}-\gamma z)^2 (\gamma-\delta \Lambda_{\rm{R}} z) \right), \nonumber \end{align} \tag{ A.12 }$$

where,

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle &C_2=\Lambda_{\rm{R}} (z (\gamma +\delta-1)-\Lambda_{\rm{R}}), \nonumber \\ &C_3=\delta z^2 \left(z^2-1\right) (\gamma +\delta -1) (\gamma z-\Lambda_{\rm{R}}) (\gamma z-\delta \Lambda_{\rm{R}}) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}}), \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle &c_5=\frac{-\delta \Lambda_{\rm{R}}^2 z^3+\Lambda_{\rm{R}} \left[\gamma +\delta +z^4 (\gamma +\delta -1)+(\gamma -1) (\delta -1) z^2-1\right)]+\gamma z (-\gamma -\delta +1)}{z^2 \left(z^2-1\right) (\gamma +\delta -1) (\gamma z-\Lambda_{\rm{R}}) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, \nonumber \\ &c_6= \frac{(1-\delta) \left(\Lambda_{\rm{R}} \left[\gamma +\delta +z^4 (\gamma\!+\!\delta\!-\!1)+(\gamma -1) (\delta -1) z^2-1\right]-\delta \Lambda_{\rm{R}}^2 z^3-\gamma z (\gamma\!+\!\delta\!-\!1)\right)}{\delta z^2 \left(z^2-1\right) (\gamma +\delta -1) (\gamma z-\Lambda_{\rm{R}}) (z (\gamma +\delta -1)-\delta \Lambda_{\rm{R}})}, \nonumber \\ &c_7=\gamma (1-\delta) \Lambda_{\rm{R}}^3 z^2 (\gamma z-\Lambda_{\rm{R}}) (z (\gamma +\delta -1)-\Lambda_{\rm{R}}) (\gamma z-\delta \Lambda_{\rm{R}}), \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \theta=& -\gamma ^2 z^6 (\gamma +\delta -1)+\delta \Lambda_{\rm{R}}^3 z \left((\delta -1) z^2+1\right)\nonumber \\ &+\Lambda_{\rm{R}}^2 \left((\delta -1)^2+z^4 (-(\gamma (\delta -1) \delta +\gamma +\delta -1))+z^2 \left(-3 \gamma \delta +\gamma -2 (\delta -1)^2\right)\right) \nonumber \\ &+\gamma \Lambda_{\rm{R}} z^3 \left(\gamma (2 \delta -1)+(\delta -1)^2+2 z^2 (\gamma +\delta -1)\right). \nonumber \end{align*}$$

A.2. Boundary vectors for Bernoulli driving

For Bernoulli driving (14) the offdiagonal left boundary vector components of the first orbital (59) are

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}_{0}} =C_1 \left(0, c_1, c_2, \chi_0 \frac{z^2 (z-2\Lambda_{\rm{L}})\Lambda_{\rm{L}}^2 (2 \Lambda_{\rm{L}}+z \mu) }{-\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}+z \mu }\right), \nonumber \end{align} \tag{ A.13 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}_{1}}=C_1 \left(0, -c_1 \left(1+\frac{1}{\mu +\nu }\right), c_3, \chi_1 \frac{z^2 (z-2 \Lambda_{\rm{L}})\Lambda_{\rm{L}}^2 (\mu +1) (2\Lambda_{\rm{L}}+z \mu)}{-\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}+z \mu }\right), \nonumber \end{align} \tag{ A.14 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{0,0}}=C_2\left(\frac{(-\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}+z \mu)^2}{\Lambda_{\rm{L}}},\Lambda_{\rm{L}} (\mu -1) \kappa, 0, 0 \right), \nonumber \end{align} \tag{ A.15 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \bra{\tilde{l}^{\prime}_{0,1}}= C_2 \left(0, 0, (z-2\Lambda_{\rm{L}})\Lambda_{\rm{L}} (\mu -1) \left(\mu z^3+\Lambda_{\rm{L}}-\Lambda_{\rm{L}} \mu \right), \frac{(-\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}+z \mu)^2}{\Lambda_{\rm{L}}} \right), \nonumber \end{align} \tag{ A.16 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle \hspace{-25mm}\bra{\tilde{l}^{\prime}_{1,0}}=C_2 \left(\frac{(2 z\Lambda_{\rm{L}}-1) (-\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}+z \mu)^2}{\Lambda_{\rm{L}}}, c_5, 0, 0 \right), \nonumber \end{align} \tag{ A.17 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\bra}[1]{{\langle #1 \vert}} \displaystyle & \bra{\tilde{l}^{\prime}_{1,1}}=C_2 \bigg(0, 0, \frac{z\Lambda_L (\mu +1) \left[\mu z^4-2\Lambda_L^2 (\mu -1) z^2+(\Lambda_L-3\Lambda_L \mu) z+4\Lambda_L^2 (\mu -1)\right]}{z-2\Lambda_L}, \nonumber \\ &\qquad\frac{(-\mu \Lambda_L+\Lambda_L+z \mu)^2}{\Lambda_L}\bigg), \nonumber \end{align} \tag{ A.18 }$$

where

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle &C_1=\frac{(\mu +1)^2}{\mu z^2 \left(z^2-1\right) (z-2\Lambda_{\rm{L}})^2 (\Lambda_{\rm{L}} (\mu -1)+z) (2\Lambda_{\rm{L}}+\mu z)^2 }, \nonumber \\ &C_2=\frac{\Lambda_{\rm{L}}^2 (\mu +1) [\Lambda_{\rm{L}} (\mu -1)-\mu ] (2 \Lambda_{\rm{L}}+\mu z -\Lambda_{\rm{L}} \mu)}{\mu z^2 \left(z^2-1\right) (2\Lambda_{\rm{L}}+\mu) \nonumber (\Lambda_{\rm{L}}+\mu z)^2}, \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_1=2 z -\Lambda_{\rm{L}} (\mu +\nu) (z-2\Lambda_{\rm{L}})^2\left(\mu z^3+2 \Lambda_{\rm{L}}\right) (z \mu -\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}}) (z-\Lambda_{\rm{L}}+\Lambda_{\rm{L}} \mu), \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_2=& -\mu ^3 z^9+\Lambda_{\rm{L}} \mu ^2 \left(-2 \mu ^2+5 \mu -1\right) z^8+2\Lambda_{\rm{L}}^2 \mu \left[6 \mu ^3-2 \mu^2+5\mu+2 (\mu +1)^2 \nu +1\right] z^7 \nonumber \\ &+\Lambda_{\rm{L}} \mu \left[2 \left(\mu ^4-10 \mu ^3-6 \mu ^2-16 \mu -8 (\mu +1)^2 \nu -5\right)\Lambda_{\rm{L}}^2+(\mu-4) \mu (\mu +1)\right] z^6 \nonumber \\ &-4\Lambda_{\rm{L}}^4\left[\mu ^5-\mu ^4-2 \mu ^2+\mu +\left(1-\mu ^2\right)^2 \nu +1\right] z^5\nonumber \\ &-\Lambda_{\rm{L}}^2\mu (\mu +1) [4 \nu +\mu (11 \mu +4 \nu -13)+10] z^5\nonumber \\ &+2\Lambda_{\rm{L}} \left(\mu ^3-\Lambda_{\rm{L}}^2 (\mu +1) \left[\mu \left\{\mu \left(\mu ^2-15 \mu -10 \nu +5\right)-10 \nu -19\right\}+2\right]\right) z^4 \nonumber \\ &+2\Lambda_{\rm{L}}^2 \left(\mu ^4-4\mu ^3+3 \mu ^2 \right)z^3\nonumber \\ &+4\Lambda_{\rm{L}}^4 (\mu +1) \left\{2 \nu +\mu (4-3 \nu +\mu [4-4 \nu+\mu (\mu +\nu -6)]-7)\right\} z^3 \nonumber \\ &+2\Lambda_{\rm{L}}^3 \left[4 \left(\mu ^2-1\right) \left(\mu ^2+\nu \mu +\nu +1\right)\Lambda_{\rm{L}}^2+\mu \left(2-5 \mu ^3+4 \mu ^2-13 \mu -2 (\mu +1)^2 \nu\right)\right] z^2 \nonumber \\ &+4\Lambda_{\rm{L}}^4 \left[(3 \mu -1) \nu (\mu +1)^2+4 \mu \left(\mu ^3+2 \mu -1\right)\right] z\nonumber \\ &+8 \Lambda_{\rm{L}}^5 \left[\nu +\mu \left(2-\mu ^3-\nu \mu ^2-\nu \mu -\mu +\nu \right)\right], \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \chi_0=&-\Lambda_{\rm{L}} \mu (\mu +1) \left(\mu ^2-2 \mu -2 \nu -1\right) z^5\nonumber \\ &+\mu (\mu -1) \left[4 \mu (\mu +1)\Lambda_{\rm{L}}^2+4 (\mu +1) \nu \Lambda_{\rm{L}}^2-\mu \right] z^4\nonumber \\ &-2\Lambda_{\rm{L}} (\mu -1)\left[(\mu +1) \left(2 \mu ^2+3 \nu \mu +\mu +\nu +1\right]\Lambda_{\rm{L}}^2-3 \mu ^2+\mu \right) z^3 \nonumber \\ &+\Lambda_{\rm{L}}^2 \left\{\mu \Big[\mu (\mu [\mu -2 \nu -17]-6 \nu +17)-2 \nu -11\Big]+2 \nu +2\right\} z^2 \nonumber \\ &-4\Lambda_{\rm{L}}^3 (\mu -1) \left(\mu ^3-4 \mu ^2-2 \nu \mu +\mu -2 \nu -2\right) z+4\Lambda_{\rm{L}}^4 (1-\mu)^2 \left(\mu ^2+\nu \mu +\nu +1\right), \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_3=&-\mu ^3 z^9+\Lambda_{\rm{L}} \mu ^2 \left(2 \mu ^2+9 \mu -1\right) z^8-2\Lambda_{\rm{L}}^2 \mu \left[6 \mu ^3+16 \mu ^2-\mu +2 (\mu +1)^2 \nu +1\right] z^7 \nonumber \\ &+\Lambda_{\rm{L}} \mu \left[2\Lambda_{\rm{L}}^2 \left(3-\mu ^4+8 \mu ^3+20 \mu ^2-2 \mu +8 (\mu +1)^2 \nu\right)-3 \mu ^3-7 \mu^2-4 \mu \right] z^6 \nonumber \\ &+\Lambda_{\rm{L}}^2 \left\{4 \nu \left[(\Lambda_{\rm{L}}-\Lambda_{\rm{L}} \mu)^2+\mu \right] (\mu +1)^2\right\}z^5\nonumber \\ &+\Lambda_{\rm{L}}^2 \left\{\mu \left[4 \left(\mu^4-2 \mu ^2+2 \mu-1\right)\Lambda_{\rm{L}}^2+\mu (17 \mu (\mu+2)+11)-6\right]\right\} z^5\nonumber \\ &+2\Lambda_{\rm{L}} \left\{\mu ^3+\Lambda_{\rm{L}}^2 (\mu +1) \left[\mu \left(\mu ^3-15 \mu ^2-11 \mu -10 (\mu +1) \nu +9\right)-2\right]\right\} z^4 \nonumber \\ &+2\Lambda_{\rm{L}}^2 \left(-\mu ^4-6 \mu ^3+3 \mu^2\right) z^3\nonumber \\ &-2\Lambda_{\rm{L}}^4 (\mu +1) \left\{2 \nu +\mu \left[-3 \nu +\mu (-4 \nu +\mu (\mu +\nu-5)-2)+4\right]-2\right\}z^3 \nonumber \\ &+2\Lambda_{\rm{L}}^3 \left[5 \mu ^3+14 \mu ^2-11 \mu -4 (\mu -1) (\mu \Lambda_{\rm{L}}+\Lambda_{\rm{L}})^2+4\right]z^2\nonumber \\ &-4\Lambda_{\rm{L}}^3\left(2\Lambda_{\rm{L}}^2 (\mu -1)-\mu \right) (\mu +1)^2 \nu z^2\nonumber \\ &-4 \Lambda_{\rm{L}}^4 \left\{\mu \Big[\nu +\mu [5 \nu +\mu (4 \mu +3 \nu +7)-7]+5\Big]-\nu -1\right\} z, \nonumber \end{align*}$$

and,

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle \chi_1=\,& (z-2\Lambda_{\rm{L}})^2 \mu \left[\Lambda_{\rm{L}} \left(\mu ^2-1\right) z^3+\mu z^2+2\Lambda_{\rm{L}} (1-\mu) z-\Lambda_{\rm{L}}^2(1-\mu)^2\right]\nonumber \\ &-2\Lambda_{\rm{L}} \nu \Big[\Lambda_{\rm{L}} \left(2 z^4-3 \Lambda_{\rm{L}} z^3-z^2+2\Lambda_{\rm{L}}^2\right) \mu ^2+\left(z^3-2 \Lambda_{\rm{L}}\right) \left(z^2-2\Lambda_{\rm{L}} z+2 \Lambda_{\rm{L}}^2\right) \mu \nonumber \\ &+\Lambda_{\rm{L}} \left(\Lambda_{\rm{L}} z^3+z^2-4 \Lambda_{\rm{L}} z+2\Lambda_{\rm{L}}^2\right)\Big], \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_4=\,&\Lambda_{\rm{L}} (\mu -1)\Big[-1-2 \nu \nonumber \\ &+\mu \left\{4 \nu z^2+\left(4 z^2-1\right) \mu -\mu[(\mu -2) \mu +2 \nu +1] z^3-2\Lambda_{\rm{L}}^2 (\mu -1)^2 (\mu +\nu)z-2 \nu -2\right\}\Big] \nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_5=\,& \mu \left(\mu ^2-2 \nu -1\right) z^3+2 \Lambda_{\rm{L}}^2 (\mu -1)^2 (\mu +\nu +1) z\nonumber \\ &-\Lambda_{\rm{L}} (\mu -1) \left\{\mu \left[4 (\mu +\nu +1) z^2-\mu -2 \nu -2\right]-2 \nu -1\right\}. \nonumber \end{align*}$$

The offdiagonal right boundary vectors components (63) are

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0,0}}=C_3\left(\frac{2 z \Lambda_{\rm{R}}}{1-2 z \Lambda_{\rm{R}}}c_6, c_6, c_7, 0 \right) \nonumber \end{align} \tag{ A.19 }$$

$$\begin{align} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0,1}}= C_3\left(0, -c_6, c_8, 0 \right) \nonumber \end{align} \tag{ A.20 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{1,0}}=C_3 \left(0, \frac{c_6 (\sigma -1)}{\sigma +1}, 0, -\frac{c_6 z \Lambda_{\rm{R}}}{\sigma +1} \right) \nonumber \end{align} \tag{ A.21 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{1,1}}=C_3\left(0, \frac{c_6}{1-2 z \Lambda_{\rm{R}}}, 0, -\frac{c_6 z \Lambda_{\rm{R}}}{\sigma +1} \right) \nonumber \end{align} \tag{ A.22 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{0}}=C_4 \left(0, c_9, c_{10}, c_9 \frac{2 \Lambda_{\rm{R}} (\rho +\sigma) (z-2 \Lambda_{\rm{R}}) (\Lambda_{\rm{R}} (\sigma-1) z+1)}{(\Lambda_{\rm{R}} (\sigma -1)+z) (2 \Lambda_{\rm{R}}+\sigma z)} \right) \nonumber \end{align} \tag{ A.23 }$$

$$\begin{align} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \renewcommand{\ket}[1]{{\vert #1 \rangle}} \displaystyle \ket{\tilde{r}_{1}}=C_4 \left(0, c_9, c_{11}, -c_9 \frac{2 \Lambda_{\rm{R}} (\rho +\sigma +1) (z-2 \Lambda_{\rm{R}}) (\Lambda_{\rm{R}} (\sigma -1) z+1)}{(\Lambda_{\rm{R}} (\sigma -1)+z) (2 \Lambda_{\rm{R}}+\sigma z)} \right) \nonumber \end{align} \tag{ A.24 }$$

where,

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle &C_3=\frac{\Lambda_{\rm{R}}^2 (\sigma +1) \left[2 \Lambda_{\rm{R}} \left(\rho \sigma +\rho +\sigma ^2+1\right)-\sigma ^2+\sigma \right] (2 \Lambda_{\rm{R}}+\sigma z)}{\sigma z^2 \left(z^2-1\right) (z-2 \Lambda_{\rm{R}}) \left[(\sigma -1) \sigma z-2 \Lambda_{\rm{R}} \left(\rho \sigma +\rho +\sigma ^2+1\right)\right]} \nonumber \\ &C_4=\frac{\sigma z^2 \left(z^2-1\right) (z-2 \Lambda_{\rm{R}}) (\Lambda_{\rm{R}} (\sigma -1)+z) \left[(\sigma -1) \sigma z-2 \Lambda_{\rm{R}} \left(\rho \sigma +\rho +\sigma ^2+1\right)\right]}{\Lambda_{\rm{R}}^3 (\sigma +1)^2 \left(2 \Lambda_{\rm{R}} \left(\rho \sigma +\rho +\sigma ^2+1\right)-\sigma ^2+\sigma \right)}\nonumber \end{align*}$$

$$\begin{align*} \renewcommand{\L}{{\rm{L}}} \newcommand{\re}{{\,{\rm{R}e}\,}} \displaystyle c_6=\,& \Lambda_{\rm{R}} (\sigma +1) z (z-2 \Lambda_{\rm{R}}) (2 \Lambda_{\rm{R}} z-1),\nonumber \\ c_7=\,&2\Lambda_{\rm{R}} (2 \rho -(\sigma -2) \sigma +1) z^4+4 \Lambda_{\rm{R}}^2 (\sigma -1) (\rho +\sigma) z^3 \nonumber \\ &-2 \Lambda_{\rm{R}} (\rho +1) (\sigma +1) z^2+[(\sigma -2) \sigma -2 \rho-1] z+2 \Lambda_{\rm{R}} [2 \rho -(\sigma -2) \sigma +1],\nonumber \\ c_8=\,& 2 \Lambda_{\rm{R}} \left(\sigma ^2-2 \rho -1\right) z^4- \nonumber \\ &4\Lambda_{\rm{R}}^2 (\sigma -1) (\rho +\sigma +1) z^3+2 \Lambda_{\rm{R}} \rho (\sigma +1) z^2+\left(1-\sigma ^2+2 \rho\right) z+2 \Lambda_{\rm{R}} \left(\sigma ^2-2 \rho -1\right),\nonumber \\ c_9=\,&\Lambda_{\rm{R}} z^2 (z-2 \Lambda_{\rm{R}}) (\Lambda_{\rm{R}} (\sigma -1)+z) (2\Lambda_{\rm{R}}+\sigma z),\nonumber \\ c_{10}=&-2 \Lambda_{\rm{R}}^2 (\sigma -1) (\sigma +1)^2-(\sigma -1) \sigma z^6+4 \Lambda_{\rm{R}} (\sigma -1) \sigma z^5+2 \Lambda_{\rm{R}}^2 (\sigma -1)z^4 \left[(\rho -1) \sigma +\rho +\sigma ^2\right]\nonumber \\ &+\Lambda_{\rm{R}} z^3 \left\{2 \rho \left[\sigma-2 \Lambda_{\rm{R}}^2 \left(\sigma ^2-1\right) +1\right]+\sigma \left[1-4 \Lambda_{\rm{R}}^2 \left(\sigma ^2-1\right)-(\sigma -2) \sigma \right]+2\right\}\nonumber \\ &-4 \Lambda_{\rm{R}}^2 z^2 \left(\rho \sigma +\rho -\sigma ^3+\sigma ^2+2\right)+4 \Lambda_{\rm{R}}^3 (\sigma -1)^2 z,\nonumber \\ c_{11}=\,& (\sigma +1)\Big(\sigma z^6+4 \Lambda_{\rm{R}}^3 (\sigma -1) z \left[z^2 (\rho +\sigma +1)-1\right] -\Lambda_{\rm{R}} z^3 \left[2 \rho +\sigma \left(1-\sigma +4 z^2\right)\right]\nonumber \\ &+2 \Lambda_{\rm{R}}^2 \left\{(\sigma +1)^2+z^4 [1-\rho \sigma+\rho -(\sigma -2) \sigma]+2 z^2 \left(\rho -\sigma^2-1\right)\right\} \Big). \nonumber \end{align*}$$