text:
\section{Introduction}

Nonequilibruim quantum dynamics has been a focus of attention in theoretical and experimental physics in the recent years. Compared to their equilibrium counterparts, nonequilibrium systems are
more challenging to describe due to the explicit time-dependence of their evolution. That generally involves coherent superposition of a large number of energy eigenstates, and not only their statistical mixture.
One key question has to do with the ways an interacting quantum system reaches \cite{deutsch1991quantum, srednicki1994chaos} (or evades \cite{moudgalya2021quantum}) thermodynamic equilibrium.
Significant progress \cite{Polkovnikov2011,nandkishore2015many,d2016quantum, Abanin2019colloquium, Bertini2021finite, Alba2021genearlized} has been achieved in understanding of the thermalization both in continuous and discrete models, especially in one-dimensional (1D) quantum systems. Some advanced numerical methods have been developed \cite{schollwock2005density, schollwock2011density, Pekker2014hilbert, cirac2021matrix}, yet they remain limited in accessible system sizes and evolution times. A number of analytical results have been obtained, especially in integrable systems  where an extensive number of conserved quantities  strongly constrain the dynamics \cite{baxter2016exactly, vsamaj2013introduction, eckle2019models}. There, the generalized hydrodynamics approach, originally introduced to study inhomogeneous quantum quenches in integrable many-body systems, has been very fruitful \cite{Bertini2021finite, Alba2021genearlized}.   Exact analytical results have  been obtained recently in systems including the Lieb-Liniger model \cite{Granet2021systematic} and the folded XXZ model \cite{pozsgay2021integrable, Zadnik2021folded1, Zadnik2021folded2}.

Typically, integrating time evolution involves discretization of time, with the small time step  ensuring that only small changes in the state of the system during a time step occur.
However, discretizing time, without assuming smallness of the time step, can also be useful for modeling physical systems.
Thus obtained quantum cellular automata (QCA) are locality-preserving maps that act in discrete space and time \cite{farrelly2020review, piroli2020quantum}, and can be regarded as a special kind of tensor networks \cite{piroli2020quantum, cirac2021matrix}. The focus of this paper are QCA that are constant-depth quantum circuits, which consist of layers of gates acting on a lattice of qubits \cite{farrelly2020review, cirac2021matrix}. Several exact results have already been obtained in such quantum circuits, including random quantum circuits \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator}, dual unitary circuits \cite{bertini2019exact, bertini2019entanglement, piroli2020exact, pavel2021correlations, claeys2021ergodic} and the Rule 54 quantum cellular automaton (QCA54)\cite{friedman2019integrable, klobas2021exact, klobas2021exact2, klobas2021ent}. The first two types are generically chaotic while QCA54 is believed to be integrable \footnote{However, it is not clear if the model can fit into a family of Yang-Baxter integrable QCA (see Ref. \cite{gombor2021integrable, prosen2021reversible} for more details).}. In QCA54, relaxation dynamics and entanglement dynamics after a quantum quench from some initial states can be computed explicitly \cite{klobas2021exact, klobas2021exact2, klobas2021ent}.

Integrable circuits can be constructed by using $\check{R}$-matrices that satisfy the Yang-Baxter equation (YBE) as the local gates  (Fig.~\ref{fig1}).  Depending on the properties of the $\check{R}$-matrix, there are three possible cases: (i) classical stochastic cellular automata (SCA) if $\check{R}$ is stochastic; (ii) unitary QCA if $\check{R}$ is unitary; (iii) nonunitary QCA if
$\check{R}$ is completely positive and trace preserving (CPTP). Most of the previous work focused on the first and the second case, while integrable nonunitary QCA are less studied.

Nonunitary quantum circuits \cite{farrelly2020review, brennen2003entanglement, piroli2020quantum}, especially unitary-projective ones \cite{li2018quantum, skinner2019meas, chan2019unitary, jian2020criticality, Ippoliti2021post, lu2021spacetime, Ippoliti2022fractal}, are important not only because they correspond to discrete-time versions of open quantum systems, but also since they can represent gate-based quantum computing and  can be used to study error-mitigation algorithms in near-term quantum devices.
One example of an integrable nonunitary QCA was recently obtained by trotterizing the Hubbard model with imaginary interaction strength \cite{sa2021integrable}. Such examples are highly nontrivial to build, which leaves a question whether there is a systematic way to construct integrable nonunitary quantum circuits.

In this paper, we construct several more examples of integrable nonunitary quantum circuits by making use of the recent developments in the systematic construction of Yang-Baxter(YB) integrable Lindblad superoperators in the continuous-time case \cite{ziolkowska2020yang, de2021constructing}. The procedure is based on finding $\check{R}$-matrices that satisfy the YBE, such that their derivatives at special points correspond to integrable Lindblad superoperators. What we demonstrate here is that thus obtained $\check{R}$-matrices not only yield continuous modes, but can also be used as the building blocks of discrete-time nonunitary quantum circtuits.
In combination with the approach of Ref.  \cite{de2021constructing} this may pave the way for a more systematic construction of nonunitary integrable quantum circuits.

The rest of the paper is organized as follows: In Section \ref{sect2}, we review the construction of integrable nonunitary quantum circuits with two-site gates and continuous integrable Lindblad superoperators. By demonstrating the CPTP property of the $\check{R}$-matrices constructed from the latter, we show that these  $\check{R}$-matrices can be used to build integrable nonunitary quantum circuits. In Section \ref{sect3}, we apply the bond-site transformation to the quantum circuits with two-site gates to construct integrable nonunitary quantum circuits with three-site gates. In Section \ref{SCA}, we briefly discuss classical SCA when the density matrix is restricted to be diagonal. In particular, we show that the Rule 150 reversible cellular automaton (RCA150) \cite{bobenko1993two} can be embedded into one of the models discussed and the Rule 105 reversible cellular automaton (RCA105) can also be embedded if the regularity condition is relaxed. In the last section, we summarize and discuss a potential way to build integrable nonunitary quantum circuits systematically, and point the way for possible future directions. Some details are relegated to the appendixes.

\section{Integrable quantum circuits with two-site gates}
\label{sect2}
\subsection{Integrability structure}
\label{int_str}
In this section, we review the integrability structure of (two-site) translationally invariant quantum circuits with two-site gates \cite{vanicat2018integrable, sa2021integrable}. It is based on the standard quantum inverse-scattering method
\cite{vsamaj2013introduction, eckle2019models, baxter2016exactly}. The central idea is to find a family of mutually commuting transfer matrices obtained from an $\check{R}$-matrix that satisfies the braid YBE. Then we can obtain a family of commuting conserved charges and even use algebraic Bethe ansatz techniques to obtain the spectrum of the transfer matrices and eigenfunctions.  For a quick introduction, see, e.g., Ref. \cite{Bertini2021finite}.

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{circuit2}
\caption{Floquet quantum circuit with period two. Each site has a Hilbert space of dimension $d$. Two-site gates $\check{R}$ that satisfy the YBE (Eq.(\ref{YBE})) act alternately.}
\label{fig1}
\end{figure}

Assume that we have an $\check{R}$-matrix, parametrized by $u$ and $v$, acting on two adjacent sites with local Hilbert space dimension $d$ each and that the $\check{R}$-matrix satisfies the braid YBE \cite{vsamaj2013introduction}:

\begin{align}
\check{R}_{23}(u, v)& \check{R}_{12}(u,w)\check{R}_{23}(v,w) =\nonumber \\
&\check{R}_{12}(v, w)\check{R}_{23}(u,w)\check{R}_{12}(u,v).
\label{YBE}
\end{align}
We let it act on each pair of sites $2j-1$ and $2j$ denoted by $(2j-1, 2j)$ with $j \in\{1, 2,..., L/2\}$ of a lattice of even length $L$ (Fig.~\ref{fig1}). We take periodic boundary conditions, namely, we identify $2L+1\sim 1$. Then we let it act on all pairs $(2j, 2j+1)$ of sites for $j \in\{1, 2,..., L/2\}$. The combination of the actions is given by
\begin{equation}
{\cal{V}}(u_1, u_2)  = \prod_{j=1}^{L/2}\check{R}_{2j, 2j+1} (u_1, u_2) \prod_{j=1}^{L/2}\check{R}_{2j-1, 2j} (u_1, u_2).
\ee
Then we repeat the action. In this way, we construct a brickwork type Floquet  quantum circuit with period 2. If we take $d =2$ for a closed system, then $\check{R}$ is 4 by 4. If we are to study open systems, it is natural to use density matrix formalism. Namely, the Floquet operator ${\cal{V}}$ acts on density matrix $\rho$. After vectorizing $\rho$, $\rho = \sum \rho_{ij} |i\rangle \langle j| \to |\rho\rangle=   \sum \rho_{ij} |i\rangle | j\rangle\rangle$, the local dimension becomes $d \to d^2$, and thus $\check{R}$ is 16 by 16.  In this section, we impose the regularity condition
\begin{equation} \check{R}(u, u) = \mathbbm{1}, \label{regcon}\ee
or equivalently, \begin{equation} R(u, u) \equiv P \check{R}(u, u) =P,\ee where $P$ is the permutation operator between two spaces defined by $P|i\rangle |j\rangle = |j\rangle |i\rangle$, which is also called the SWAP gate. It describes the trivial exchange process of the scattering of two particles \cite{vsamaj2013introduction}.  This condition, together with the YBE, implies \cite{gombor2021integrable} the inversion property,
\begin{equation} \check{R}_{12}(u, v) \check{R}_{12}(v, u) =\mathbbm{1}.\label{inv}\ee

To show that the circuit is YB integrable, we follow Ref.~\cite{sa2021integrable} and define a staggered transfer matrix:
\begin{align}
T(u) & =  \nonumber \\
&\text{Tr}_{a}[R_{a L} (u, u_1 )R_{a L-1} (u, u_2 )... R_{a 2} (u, u_1 ) R_{a 1} (u, u_2 )]
\end{align}
where $a$ can be interpreted as an auxiliary space to take trace over. Then the YBE leads to
$[T(u), T(v)]=0$ for any $u$ and $v$.

The regularity property and the inversion property then imply
\begin{align}
 &T(u_1)   \nonumber \\
=&P_{12} P_{13}...P_{1L}   \check{R}_{12} (u_1, u_2 )\check{R}_{34} (u_1, u_2 )... \check{R}_{L-1 L} (u_1, u_2 ) \nonumber \\
 \equiv & U \check{R}_{12} (u_1, u_2 )\check{R}_{34} (u_1, u_2 )... \check{R}_{L-1 L} (u_1, u_2 ),
\end{align}
and
\begin{align}
& T(u_2) \nonumber \\
=& P_{12} P_{13}...P_{1L} \check{R}_{23}^{-1} (u_1, u_2 )\check{R}_{45}^{-1} (u_1, u_2 )... \check{R}^{-1}_{L 1} (u_1, u_2 ) \nonumber \\
 \equiv & U  \check{R}_{23}^{-1} (u_1, u_2 )\check{R}_{45}^{-1} (u_1, u_2 )... \check{R}^{-1}_{L 1} (u_1, u_2 ),
\end{align}
where $U =P_{12} P_{13}...P_{1L}$ is a one-site shift. We have used the identity ${\cal{O}}_{ab} P_{ac} = P_{ac} {\cal{O}}_{cb}$ (for $b\neq c$) and $\text{Tr}_a P_{a1} = \mathbbm{1}$. Thus,

\begin{equation}
{\cal{V}} (u_1, u_2)  = T(u_2)^{-1} T(u_1)
\ee
and it commutes with $T(u)$ for any $u$. The commuting charges are given by
\begin{equation}
Q_n^+ =\frac{d^n}{d u^n}  \ln T(u) |_{u =u_1}, \quad Q_n^- =\frac{d^n}{d u^n}  \ln T(u) |_{u =u_2}.
\ee

\subsection{Integrable Lindblad Superoperators}\label{ILS}
Quantum circuits that we discussed may be regarded as coming from trotterization of  corresponding lattice models. If an open continuum system is in contact with a Markovian environment whose dynamics is much faster than  the internal one, the equation of motion is approximately given by the Lindbladian master equation
\begin{equation}
\dot{\rho} = {\cal{L}} \rho = i[\rho, H] + \sum_a l_a \rho l_a^{\dagger} - \frac{1}{2}\{ l_a^{\dagger} l_a, \rho\},
\label{lme}
\ee
where $\rho$ is the density operator, ${\cal{L}}$ is the Lindbladian and $l_a$ are jump operators. For a lattice system,  ${\cal{L}} = \sum_i {\cal{L}}_{i, i+1}$.

A recent study \cite{de2021constructing} constructed many new examples of integrable Lindblad superoperators following Ref.~\cite{de2019classifying, de2020classifying}.  We briefly summarize their approach. In terms of the vectorized density matrix, $\rho \in {\cal{H}}^{(1)} \otimes  {\cal{H}}^{(2)} $,  the Lindbladian density for a system with one family of jump operators has the explicit form
\begin{align}
{\cal{L}}_{i, i+1}(u) & = -i h_{i, i+1}^{(1)}(u) + i h_{i, i+1}^{(2)*}(u) +l_{i, i+1}^{(1)}(u) l_{i, i+1}^{(2)*}(u)  \nonumber \\
& - \frac{1}{2} l_{i, i+1}^{(1)\dagger}(u) l_{i, i+1}^{(1)}(u) - \frac{1}{2} l_{i, i+1}^{(2)\text{T}}(u) l_{i, i+1}^{(2)*}(u).
\label{lin}
\end{align}
Here, ${\cal{H}}^{(1)}$ is that of two qubits and ${\cal{H}}^{(2)}$ is its canonical dual space; $h^{(s)}$ and $l^{(s)}$ act on the $s$-th Hilbert space. $^*$ denotes the complex conjugate and $^\text{T}$ transpose. The dependence on the spectral parameter $u$ will be used in the following.

To construct integrable ${\cal{L}}$, we use ansatze for ${\cal{L}}$ of the form in Eq.(\ref{lin}) and its associated $\check{R}$ satisfying
\begin{equation}
\check{R}(u, u) =\mathbbm{1}, \quad   {\cal{L}}_{i, i+1}(u) = \partial_v \check{R}_{i, i+1}(v, u)|_{v\to u}.
\label{initial}
\ee

Let $Q_n$ be the $n$-site conserved charge, in particular, $Q_2 = {\cal{L}}$. Higher charges are generated by
\begin{equation}
Q_{n+1}(u) = [ B(u), Q_n(u)], \quad n >1,
\label{q2}
\ee
where the boost operator $B(u):=  \partial_u + \sum_{i = -\infty}^{\infty} i {\cal{L}}_{i, i+1}(u)$. This expression is defined for infinite chains but reduces to finite chains consistently. Then impose the condition
\begin{equation} [Q_2(u), Q_3(u)] =0\ee on the ansatz and solve the following equations obtained from the YBE:
\begin{subequations}
\begin{gather}
[R_{13} R_{23}, {\cal{L}}_{12}(u)] = \dot{R}_{13} R_{23} - R_{13} \dot{R}_{23};\\
[R_{13} R_{12}, {\cal{L}}_{23}(v)] = R_{13} R'_{12} - R'_{13} R_{12}.
\end{gather}
\label{suther}
\end{subequations}
In the first equation $R_{ij} = R_{ij}(u, u)$ and $\dot{R}$ is the derivative taken with respect to the first varaible at $v=u$. In the second equation, $R_{ij} = R_{ij}(v, v)$ and the derivative is taken with respect to the second variable at $u =v$. Combining Eqs.(\ref{initial}), (\ref{q2}) and (\ref{suther}), we can obtain $\cal{L}$ and $\check{R}$. In the end, we can check if the $\check{R}$-matrix obtained this way satisfies the YBE (Eq.(\ref{YBE})).

Note that there are some degrees of freedom in $h$ and $l$ that keep ${\cal{L}}_{i, i+1}(u)$ invariant \cite{de2021constructing}, thus yielding the same $\check{R}$. However, different ${\cal{L}}_{i, i+1}(u)$ may give rise to $\check{R}$-matrices that are still equivalent via local basis transformations, reparameterizations, normalizations, discrete transformations, and twists. See Ref.\cite{de2019classifying, de2020classifying} for more information. In order to qualify as integrable non-unitary (open) quantum circuits,  the $\check{R}$-matrix has to be a quantum channel or, equivalently, satisfy CPTP \cite{preskill1998lecture}.

\subsection{CPTP}
To show that an $\check{R}$-matrix is CPTP, we can use the Choi-Jamiolkowski representation \cite{sa2021integrable} of $\check{R}$ and define the Choi matrix $D_{\beta \delta \zeta \theta}^{\alpha \gamma\varepsilon \eta} \equiv \check{R}_{\gamma \delta \eta \theta}^{\alpha \beta \varepsilon \zeta}$. Then, $\check{R}$ is trace-preserving if and only if  $D^{\alpha \gamma \varepsilon\eta}_{\alpha \delta \varepsilon \theta} = \check{R}^{\alpha\alpha \varepsilon\varepsilon}_{\gamma \delta \eta \theta} = \delta_{\alpha}^{\gamma} \delta_{\theta}^{\eta}$; $\check{R}$ is completely positive if and only if $D$ is non-negative. Additionally,  $\check{R}$ is unital if and only if $D^{\alpha \gamma \varepsilon\eta}_{\beta \gamma \zeta \eta} = \check{R}^{\alpha\beta \varepsilon\zeta }_{\gamma \gamma \eta \eta} = \delta^{\alpha}_{\beta} \delta^{\varepsilon}_{\zeta}$. Unitality is not necessary for $\check{R}$ to be CPTP.

It is easy to check explicitly that all $\check{R}$-matrices of all five cases in Ref.~\cite{de2021constructing} are CPTP. Additionally,  \begin{equation} \text{tr}(D (u))= \text{tr}(D(0)) =4.
\ee
Hence, it has at most rank 2 and the two non-vanishing eigenvalues satisfy
\begin{equation}
\lambda_1 +\lambda_2 =4,
\ee
and $\lambda_1, \lambda_2 \ge 0$ for $u\ge 0$. Since $\check{R}$ is CPTP, we have a Kraus representation \cite{preskill1998lecture} of  $\phi_{\beta \delta \zeta \theta}^{\alpha \gamma\varepsilon \eta} \equiv \check{R}_{\beta \zeta \delta \theta}^{\alpha \varepsilon \gamma \eta}$:
\begin{equation}
\phi  =K_0\otimes K_0^* + K_1\otimes K_1^* ,
\ee
which satisfies
\begin{equation}
K_0^{\dagger}K_0 +K_1^{\dagger}K_1 = \mathbbm{1}.
\ee
It acts on $\rho$ as
\begin{equation}
\phi[\rho] = K_0\rho K_0^{\dagger} + K_1 \rho K_1^{\dagger}.
\ee
$\check{R}$ is unital if
\begin{equation}
K_0 K_0^{\dagger} +K_1  K_1^{\dagger}= \mathbbm{1}.
\label{unital}
\ee

Thus, we can use these $\check{R}$-matrices to construct integrable nonunitary quantum circuits.
Note that the CPTP property of $\check{R}$ holds for all $u \ge 0$. However, when $u\to 0$, the quantum circuits reduce to the continuous systems described by the Lindbladian master equation (Eq.(\ref{lme})). Indeed, in a closed system with $d=2$, the $u\to 0$ limit of the quantum circuit described in Subsection \ref{int_str} reduces to $\lim_{N\to \infty} {\cal{V}}(u/N)^N  \to \exp(i u H)$, where $H$ can be interpreted as the hermitian Hamiltonian \cite{vanicat2018integrable}. Similarly, in the open system, when we take $u\to 0$, $ \lim_{N\to \infty} {\cal{V}}(u/N)^N \to \exp(u{\cal{L}})$.
The quantum circuit reduces to a system with continuous `time' $u$ \cite{vanicat2018integrable}.

In general, if there are more than one family of jump operators, the expansion of the Kraus operators around $u =0$ is  \footnote{Note that if $\check{R}$ has a difference form, i.e., $\check{R}(u, v) = \check{R}(u-v)$, we simply write $\check{R}(u) = \check{R}(u, 0)$. Similarly for $K_a$ and other operators. In this case, ${\cal{L}}$ is independent on $u$.}
\begin{subequations}
    \begin{align}
K_0(u) & = I + (-i H - \frac{1}{2} \sum_{a >0} l_a^{\dagger} l_a ) u + O(u^2), \\
K_a (u) & = \sqrt{u} l_a  + O(u).
\end{align}
\label{kraus}
\end{subequations}

In the following we show that the five $\check{R}(u)$-matrices in  Ref.\cite{de2021constructing} constructed using integrable Lindblad superoperators are CPTP. Indeed, it would be interesting to see if all $\check{R}$-matrices constructed in this way are always CPTP; so far we have not been able to demonstrate that.
In principle, we can also run the argument in reverse: If we have a YBE solution $\check{R}(u)$ that has the CPTP property, it can be used to obtain a corresponding integrable Lindblad superoperator.

\subsection{Examples}
We list all five cases in
Ref.~\cite{de2021constructing}. We can classify them into two categories. The first category includes Model B1 and Model B2 where two limiting cases ($u\to 0$ and $u\to \infty$) are both unitary. The rest (Model A1, Model A2, and Model B3) do not have this property.
\subsubsection{Model B1}
The (16 by 16) $R$-matrix we use is obtained from that given in Ref.~\cite{de2021constructing} by rescaling it by a factor $e^u/(u+1)$. $h =0$ and
\begin{equation}
l =
\begin{pmatrix} \tau & 0 & 0 & 0\\
0 & 0& 1 & 0 \\
0 &1 &0  & 0\\
0 & 0 & 0 & \kappa
\end{pmatrix},
\ee
where $\tau = \pm 1$ and $\kappa =\pm 1$. For $\tau =\kappa =1$,
\begin{equation}
\check{R}(u) = \frac{\mathbbm{1} +u P}{u+1},
\ee
where $P$ is the permutation matrix.
This is the $R$-matrix of the $SU(4)$ invariant chain if $u$ is taken to be purely imaginary. It is easy to see that $\check{R} = P R$ satisfies the inversion property Eq.(\ref{inv}).

In this case, the $\check{R}$-matrix is both CPTP and unital. The non-trivial eigenvalues of the Choi matrix $D$ are
\begin{equation}
\lambda_{\pm} = 2 \pm 2 \frac{\sqrt{u^2 -(\frac{3-\tau \kappa}{2} )u+1}}{u+1}.
\ee

The Kraus operators are
\begin{equation} K_0(u) = \sqrt{\frac{1}{u+1}} \mathbbm{1}_{4\times 4}, K_1(u) = \sqrt{\frac{u}{u+1}}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& \tau & 0 \\
0 &\tau &0  & 0\\
0 & 0 & 0 & \kappa\tau
\end{pmatrix}.
\ee

The action on $\rho$ is
\begin{equation}
\phi[\rho] = \frac{\rho + u K_1 \rho K_1 }{u+1}.
\label{Kraus}
\ee
It is easy to check that  $K^{\dagger}_0 K_0 + K_1^{\dagger} K_1 =\mathbbm{1}$.  When $u =0$, $K_0(0) = \mathbbm{1}$ and $K_1(0) =0$. When $u\to \infty$,
\begin{equation} K_0(\infty) = 0, \quad K_1(\infty) =
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& \tau & 0 \\
0 &\tau &0  & 0\\
0 & 0 & 0 & \kappa\tau
\end{pmatrix}.
\label{k0k1}
\ee
Thus, the trivial unitary quantum circuit $\check{R} =\mathbbm{1}$ is connected to the unitary circuit $U = K_1(\infty)$ by tuning $u$ from 0 to $\infty$.  For $\tau \kappa =1$, $K_1(\infty)$ corresponds to the XXZ model with $\Delta =\pm 1$. In this case, $K_1(\infty)$ will not generate entanglement for any (pure) product states, which is not true for the other two cases. For $\tau \kappa =- 1$, it corresponds to the XX model with a specific $h_z$. For all four cases, $K_1(\infty)$ is dual-unitary, and it belongs to the noninteracting case described in Ref. \cite{bertini2019exact}.

That $\check{R}$ is unital implies that the maximally mixed state is a steady state. It is the quantum version of the bi-stochastic classical map \cite{pavel2021correlations}. The conserved charges can be obtained as in Ref.\cite{vanicat2018integrable}.

\subsubsection{Model B2}
This case includes the XX model with dephasing noise (which corresponds to the Hubbard model with imaginary coupling) discussed in Ref.~\cite{sa2021integrable} as the special case ($\phi =\pi/2$ in the following discussion). $h$ and $l$, respectively, are given below:
\begin{subequations}
\begin{equation}
h =2  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\begin{align}
&\frac{l}{\eta} = \nonumber \\
& \begin{pmatrix} \cosh(2u) & 0 & 0 & 0\\
0 & 1& \text{i} e^{\text{i} \phi} \sinh(2u) & 0 \\
0 & -\text{i} e^{-\text{i} \phi} \sinh(2u)& -1  & 0\\
0 & 0 & 0 & -\cosh(2u)
\end{pmatrix},
\end{align}
\end{subequations}
where \begin{equation} \eta (u) = \sqrt[4]{\frac{w^2}{ w^2\sinh^2 (2u) +1}}.\ee

Define
\begin{equation}
\check{G}(u) =\frac{1}{a}  \begin{pmatrix} a & 0 & 0& 0\\ 0 & c &-\text{i} b e^{\text{i} \phi} & 0 \\  0 & -\text{i} b e^{-\text{i} \phi} & c& 0 \\ 0 & 0 & 0& a\end{pmatrix},
\ee
where \begin{equation} a = \cosh (u), \quad b = \sinh(u), \quad c=1.\ee
Also, take the matrix $e^{\beta}_{\alpha}$ whose matrix element is $(e^{\beta}_{\alpha})^{\beta'}_{\alpha'} =\delta_{\beta, \beta'}\delta_{\alpha, \alpha'}$, and
\begin{subequations}
\begin{gather}
\check{r}_{\uparrow}(\lambda) = \check{G}^{\alpha \gamma}_{\beta \delta}(\lambda) e_{\alpha}^{\beta} \otimes \mathbbm{1}_2 \otimes e_{\gamma}^{\delta} \otimes \mathbbm{1}_2, \\
\check{r}_{\downarrow}(\lambda) = \check{G}^{\alpha \gamma}_{\beta \delta}(\lambda)\mathbbm{1}_2\otimes e_{\alpha}^{\beta} \otimes \mathbbm{1}_2 \otimes e_{\gamma}^{\delta},
\end{gather}
\end{subequations}
Define $\check{r} = \check{r}_{\uparrow} \check{r}^*_{\downarrow}$ and then the $\check{R}$-matrix can take the form
\begin{equation}
\check{R}(u, v) = \beta \check{r}(u-v) + \alpha \check{r}(u+v) (\sigma^z \otimes \sigma^z \otimes \mathbbm{1}_4),
\ee
where
\begin{equation}
\beta =
  \frac{\cosh(u - v) \cosh(
     \xi - \eta)}{(\cosh (u - v) \cosh( \xi - \eta) + \cosh(u + v) \sinh( \xi - \eta))},
\ee
$\alpha = 1-\beta$, and
\begin{equation}
\frac{\sinh(2\xi)}{\sinh(2u)} =\frac{\sinh(2\eta)}{\sinh(2v)} = w, \quad u\ge v.
\ee
The inversion condition, Eq.(\ref{inv}), is satisfied. $\check{R}(u, v)$ is not of difference form  and it is unital and CPTP. In particular, the nonvanishing eigenvalues of its Choi matrix are $4\alpha$ and $4\beta$. The Kraus operators are given by
\begin{equation}
K_0 = \sqrt{\beta}   \check{G}(u-v),\quad  K_1 = \sqrt{\alpha}   \check{G}(u+v) (\sigma^z \otimes \mathbbm{1}).
\ee
When $u=v$, $\xi =\eta$, $\beta =1$ and $\alpha=0$. $K_0(0)=\mathbbm{1}$ and $K_1(0) =0$. When $u \to \infty$, $\alpha =1$ and $\beta=0$. Then
\begin{equation} K_0(\infty) = 0, \quad K_1(\infty) =
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& -i e^{i \phi} & 0 \\
0 &i e^{-i \phi} &0  & 0\\
0 & 0 & 0 & -1
\end{pmatrix}.
\ee
Note that when $\phi =\pi/2$ this is same as Eq.(\ref{k0k1}) for $\tau=-1$ and $\kappa =-1$. Both limiting cases are unitary. The integrability of the circuit is also demonstrated numerically in Ref.\cite{sa2021integrable} by examining the spectral statistics.

Note that both Model B1 and Model B2 have the following form
$\phi[\rho] = p U_0\rho U_10^{\dagger} + (1-p) U_1\rho U_1^{\dagger}$ with $0\le p \le 1$.  More generally, $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$ has the form of a convex combination where $U_i$ is unitary. It resembles random unitary quantum circuits if we take $U_i$ to be a Haar-random unitary \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator} or  Clifford-random unitary \cite{li2018quantum, Ippoliti2021post}.

\subsubsection{Model A1}
In this model, $h$ and $l$ are given by
\begin{equation}
h =\frac{1}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},\quad
l=\\
 \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& 0& 0 \\
0 & 1& -i e^{-\text{i} \phi}  & 0\\
0 & 0 & 0 &0
\end{pmatrix}.
\ee
The entries of the $R$-matrix are given in Ref. \cite{de2021constructing}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital in this case since Eq.(\ref{unital}) is not satisfied. The two nonvanishing eigenvalues of the Choi matrix are $
\lambda_{\pm} =2 \pm \sqrt{1 +2 e^{-2 u} + e^{-3 u}}$. Thus
$\lambda_{\pm} \ge 0$ when $u \ge 0$.
The Kraus operators are given by
\begin{subequations}
    \begin{align}
 K_0 = e^{-u/2}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & e^{-u/2} & 0 & 0 \\
0 &0 &e^{u/2}  & 0\\
0 & 0 & 0 & 1
\end{pmatrix}, \\
 K_1 =  \sqrt{1-e^{-u}}
\begin{pmatrix} 1& 0 & 0 & 0\\
0 & e^{-u/2} & 0 & 0 \\
0 &-\text{i} e^{ \text{i} \phi} & 0 & 0\\
0 & 0 & 0 & 1
\end{pmatrix}.
\end{align}
\end{subequations}
The phases of $K_0$ and $K_1$ are irrelevant.

\subsubsection{Model A2}
Here, $h$ and $l$ are given by
\begin{equation}
h =\frac{\tau}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& i  & 0 \\
0 &-i &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},\quad
l=\\
 \begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& 0& 0 \\
0 & \tau & 1  & 0\\
1 & 0 & 0 &0
\end{pmatrix},
\ee
where $\tau =\pm 1$.
The entries of the $R$-matrix of the Choi matrix are given in Ref. \cite{de2021constructing}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital in this case since Eq.(\ref{unital}) is not satisfied.
The two nonvanishing eigenvalues are
$\lambda_{\pm} = 2 \pm 2 e^{-3u/2}$ and $\lambda_{\pm} \ge 0$ for $u \ge 0$. The Kraus operators are given by
\begin{subequations}
    \begin{align}
    K_0= e^{-u}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & e^{u/2} & 0 & 0 \\
0 &-\tau (e^u-1)  &e^{u/2}  & 0\\
0 & 0 & 0 & e^{u}
\end{pmatrix}, \\ K_1 = \sqrt{e^u -1} e^{-u}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
0 &\tau &e^{u/2}  & 0\\
e^{u/2} & 0 & 0 & 0
\end{pmatrix}.
\end{align}
\end{subequations}
 $\check{R}$ is not unital in this case since Eq.(\ref{unital}) is not satisfied.

\subsubsection{Model B3}
$h$ and $l$ are given by
\begin{subequations}
 \begin{equation}
h  =\frac{1}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\begin{equation}
l =\sqrt{\frac{\gamma}{2}}
 \begin{pmatrix} \gamma & 0 & 0 & 0\\
0 & 1& i (\gamma-1) e^{i \phi}& 0 \\
0 & -i (\gamma+1) e^{-i \phi}& -1  & 0\\
0 & 0 & 0 &\gamma
\end{pmatrix}.
\ee
\end{subequations}
$\gamma\ge 0$.
The entries of the $R$-matrix of the Choi matrix are given in Ref. \cite{de2021constructing}   \footnote{There are typos for the $R$-matrix of Model B3 in Ref. \cite{de2021constructing} : $R_4^{10} = e^{\text{i} \phi} 2\text{i} (\gamma+1)^2 \sinh(\alpha)/\zeta$ and $R_4^7 = -e^{-\text{i} \phi} 2\text{i} (\gamma+1)^2 \sinh(\alpha)/\zeta$.}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital for generic $\gamma$ in this case since Eq.(\ref{unital}) is not satisfied.  The Kraus operators are given by

\begin{subequations}
    \begin{align}
   &  K_0 = \nonumber \\
     & \frac{1}{\eta}
\begin{pmatrix} \eta  & 0 & 0 & 0\\
0 & 1 & -\text{i} (1-\gamma)\sinh (\alpha)e^{i \phi} & 0 \\
0 & -\text{i} (1+\gamma)\sinh (\alpha) e^{-i \phi} & 1 & 0\\
0 & 0 & 0 &  \eta
\end{pmatrix},
\end{align}
\begin{equation}  K_1 = \frac{\xi}{\eta}
\begin{pmatrix} 0 & 0 & 0 & 0\\
0 & -e^{-\alpha} (1-\gamma) & i(1-\gamma) e^{i \phi}& 0 \\
0 & i(1+\gamma)e^{-i \phi} &  e^{\alpha} (1+\gamma) & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\end{subequations}
where $\alpha  = (\gamma^2+1)u/2$, $\eta  = \gamma \sinh(\alpha) +\cosh(\alpha)$, and $\xi= \sqrt{\gamma \sinh(\alpha) /[(1+\gamma^2) \cosh(\alpha) +2\gamma \sinh(\alpha)] }$.

\section{Integrable quantum circuits with three-site gates}
\label{sect3}

In the previous section we showed how two-site $\check{R}$ matrices that appear in the integrable Lindbladian construction can be repurposed to build non nonunitary quantum circuits. In this section, we extend this construction to integrable three-site nonunitary gates and circtuits by taking advantage of bond-site transformations.

\subsection{Bond-site transformation}
\label{sec_bst}
A bond-site transformation (BST) maps a relative orientation of a pair of spin to a single  spin value.  We discuss BST following Ref.~\cite{gombor2021integrable}. For a 1D spin chain, we choose the computational basis $\mid s_1, s_2, ..., s_L\rangle$ where $s_i =\pm 1$. For each pair of neighboring spins, define a state $s'_i = s_i s_{i+1}$ on the bond. More concretely, if two neighboring sites are in the same state, then $s' =1$; if two neighbor sites are in different states,then  $s' =-1$. Thus, the state $s'=-1$ can be interpreted as the existence of a domain wall. Note that this is a two-to-one correspondence: states related by a global spin reflection are mapped to the same state on the bond. If open boundary conditions are imposed on the original spin chain and one of the boundary states is fixed, the states on the bonds after the BST are in the one-to-one correspondence to the rest of the states. If periodic boundary conditions are taken, the mapping is two-to-one, and the bond states are constrained by $s'_1 s'_2\cdots s'_L =1$. Note that, similar to the Jordan-Wigner transformation \cite{fradkin2013field}, a BST is nonlocal because a domain wall is nonlocal in the original basis. Moreover, an entangled state can be mapped to an unentangled one. For example, the Greenberger–Horne–Zeilinger state $(\mid\uparrow\rangle^{\otimes L} +\mid\downarrow\rangle^{\otimes L})/\sqrt{2}$ is mapped to a product state.

On the operator level, an $(l+1)$-site operator may be mapped to an $l$-site operator after a BST. In particular, three-site operators that commute with the global spin reflection ($\Sigma^x_1\otimes \Sigma^x_2 \otimes \Sigma^x_3$) are in one-to-one correspondence with a two-site operators that commute with $\sigma^z \otimes \sigma^z$.

\begin{figure}[H]
\centering
\includegraphics[width=0.2\textwidth]{b.s.t.1}
\label{figbts}
\end{figure}

They are listed below:
\begin{align}
\sigma^x & \otimes \sigma^x  \leftrightarrow \mathbbm{1} \ \otimes \Sigma^x \otimes\ \mathbbm{1}, \quad
\sigma^y \otimes \sigma^y  \leftrightarrow - \Sigma^z \otimes \Sigma^x \otimes \Sigma^z,\nonumber  \\
\sigma^z& \otimes \sigma^z  \leftrightarrow \Sigma^z \otimes\ \mathbbm{1} \ \otimes \Sigma^z, \quad
\sigma^y   \otimes \sigma^x  \leftrightarrow \Sigma^z \otimes \Sigma^y \otimes\ \mathbbm{1},\nonumber \\
\sigma^x &\otimes \sigma^y  \leftrightarrow\ \mathbbm{1} \  \otimes \Sigma^y \otimes \Sigma^z, \quad
\mathbbm{1}   \otimes \sigma^z  \leftrightarrow -\mathbbm{1}\ \otimes \Sigma^z \otimes \Sigma^z,\nonumber \\
\sigma^z &\otimes \mathbbm{1} \leftrightarrow -\Sigma^z \otimes \Sigma^z \otimes\ \mathbbm{1}.
\end{align}
Consequently, general operators of the following form
\begin{align}
\check{R} &= h_1 \mathbbm{1} \otimes \mathbbm{1} + \frac{h_2}{2} (\sigma^z \otimes \mathbbm{1} - \mathbbm{1} \otimes \sigma^z) + h_3 \sigma^+ \otimes \sigma^- \nonumber \\
&+ h_4 \sigma^- \otimes \sigma^+ + \frac{h_5}{2} (\sigma^z \otimes \mathbbm{1} + \mathbbm{1} \otimes \sigma^z) + h_6 \sigma^z \otimes \sigma^z \nonumber \\
&+ h_7 \sigma^-\otimes \sigma^- + h_8 \sigma^+ \otimes \sigma^+.
\end{align}
can all be mapped to three-site operators.
Note that this is the general form of the ansatz used in Ref.~\cite{de2020classifying}. After a BST, the corresponding three-site operator $\check{G}$ satisfies
\begin{equation}
[\check{G}_{123}, \check{G}_{345}]=0,
\label{comm}
\ee
since the operator acting on the common site (3) is either $\mathbbm{1}$ or $\Sigma^z$.

Since we want to consider open quantum systems, we need to generalize the transformation to mixed density matrices.
The vectorized density matrix, $|\rho\rangle$, is generated by $|s'_1,...\rangle |\tilde{s}'_1,...\rangle\rangle $. In particular, the basis states upon which the two-site gate $\check{R}$ discussed in the previous section acts are generated by $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$.  Since we constructed $\check{R}$ using the Lindbladian of the form in Eq.(\ref{lin}), if both $h$ and $l$ commute with $\sigma^z\otimes \sigma^z$, the action of $\check{R}$ on $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ may preserve the product of $s'_1 s'_2$ and of $\tilde{s}'_1 \tilde{s}'_2$, separately. If this is the case, $|s'_1, s'_2\rangle$ and $ |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ can separately be mapped to $|s_1, s_2, s_3\rangle$ and $|\tilde{s}_1, \tilde{s}_2, \tilde{s}_3\rangle\rangle$. There are four states corresponding to each $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ and the (64-dim) Hilbert space separates into four sectors. Moreover, $\check{R}$ can be regarded as a linear operator within each of the sectors. The corresponding three-site gate $\check{G}$ then also preserves each sector. Thus, $\check{G}$ is equal to the direct sum of four copies of $\check{R}$ when acting on the (64-dim) vectorized density operator, which, by linearity, implies that
\begin{align}
\check{R}_{23}(u_1, u_2)&\check{R}_{12}(u_1, u_3) \check{R}_{23}(u_2, u_3) \nonumber \\
&=\check{R}_{12}(u_2, u_3)\check{R}_{23}(u_1, u_3)\check{R}_{12}(u_2, u_3).
\label{RYBE}
\end{align}
leads to
\begin{align}
\check{G}_{234}(u_1, u_2)& \check{G}_{123}(u_1, u_3) \check{G}_{234} (u_2, u_3) \nonumber \\
& = \check{G}_{123}(u_2, u_3) \check{G}_{234}(u_1, u_3) \check{G}_{123} (u_1, u_2).
\label{GYBE}
\end{align}
Also, if two-site (Kraus) operators $K_i$ satisfy
$
K_0^{\dagger}K_0 +K_1^{\dagger}K_1 = \mathbbm{1}
$, the corresponding three-site operators $\tilde{K}_i = K_i\otimes K_i$ satisfy
$
\tilde{K}_0^{\dagger}\tilde{K}_0 +\tilde{K}_1^{\dagger}\tilde{K}_1 = \mathbbm{1}.
$

\begin{figure}[b]
\centering
\includegraphics[width=0.4\textwidth]{circuit3}
\caption{Floquet quantum circuit with staggered special three-site gates. A blue disk represents a general operator while a black dot can be $\mathbbm{1}$ or $\Sigma^z$.}
\label{circuit3}
\end{figure}

\subsection{Integrability structure}
Since the BST is nonlocal, it is  not obvious whether all  local charges of a circuit with two-site gates are mapped to local charges of the corresponding circuit with three-site gates. To show that this is indeed the case, and thus the latter is integrable, we follow Ref. \cite{gombor2021integrable}.

We have already seen that if the two-site gate $\check{R}$ satisfies the YBE (Eq.(\ref{RYBE})) the corresponding three-site gate $\check{G}$ satisfies Eq.(\ref{comm}), Eq.(\ref{GYBE}), $ \check{G}_{123}(v, v)= \mathbbm{1}$, and
\begin{equation} \check{G}_{123}(u, v) \check{G}_{123}(v, u) = \mathbbm{1}.\ee
Define $\check{L}_{123} (u) =\check{G}_{123}(u, 0)$ and
\begin{equation}
\check{R}_{12, 34}(u,v) = \check{L}_{234}(u) \check{G}_{123}(u,v) \check{L}_{234}(v)^{-1}.
\label{rlgl}
\ee
Thus defined four-site $\check{R}$ then satisfies
\begin{align}
&\check{R}_{34, 56}(u_1, u_2) \check{R}_{12, 34}(u_1, u_3) \check{R}_{34, 56} (u_2, u_3) \nonumber \\
& = \check{R}_{12, 34}(u_2, u_3) \check{R}_{34, 56}(u_1, u_3) \check{R}_{12, 34} (u_1, u_2),
\end{align}
and we can construct a transfer matrix
\begin{align}
T(u) &= \text{Tr}_{a,b} [R_{(1,2), (a,b)}(v, u)  R_{(3,4), (a,b)}(v, u)... \nonumber \\
& R_{(L-3,L-2), (a,b)}(v, u) R_{(L-1,L), (a,b)}(v, u)],
\end{align}
where $(a, b)$ denote the ``auxiliary space" and $R_{(a,b), (c,d)} \equiv P_{(a,b), (c,d)}\check{R}_{(a,b), (c,d)}$ ($P$ is the permutation operator).
Regarding each pair of indices as a single index, we can follow the derivation Section \ref{int_str} to show that $ T(v)  = U^{-2}$, where $U$ is the one-site translation operator, and
\begin{equation} [T(u), T(u')] =0.\ee
Moreover, by Eq.(\ref{rlgl}),  $\check{R}_{12, 34}(u,0) = \check{L}_{234}(u) \check{L}_{123}(u)$. Then
\begin{align}
 T(0) &=  \nonumber \\
&\text{Tr}_{a,b} [L_{1,2, b}L_{1,2, a}  L_{3,4, b} L_{3,4, a}\cdots   L_{L-1,L, b}L_{L-1,L, a}],
\end{align}
Define
\begin{equation} {\cal{V}}(v, u) =  T(v)^{-1} T(u),\ee
then the Floquet operator (Fig.\ref{circuit3}) corresponds to
\begin{equation} {\cal{V}}(v, 0) = \prod_{j=1}^{L/2}\check{L}_{2j, 2j+1, 2j+2}(v)\prod_{j=1}^{L/2}\check{L}_{2j-1, 2j, 2j+1}(v).\ee
The operator $ {\cal{V}}(v, u) $  leads to local conserved charges, the first being
\begin{equation}
Q_4  =\frac{d }{d u }  {\cal{V}}(v,u) |_{u =v}.
\ee
It is a sum of four-site operators. Higher derivative yield higher conserved charges. Therefore, we conclude that, if the BST is applicable to an integrable Floquet quantum circuit with two-site gates as discussed above, the corresponding circuit with three-site gates is also integrable. In particular, all five models discussed in the last section can be mapped to integrable circuits with nonunitary three-site gates.

\subsection{Example}
\label{example}
Let us discuss the simplest example:  Model B1 with $\tau \kappa =1$. Both Kraus operators of Model B1 are special cases of a general $\check{R}$-matrix of the zero-field eight-vertex (ZFEV) model:
\begin{equation}
\check{R} =  \begin{pmatrix} a & 0 & 0& d\\ 0 & c & b& 0 \\  0 & b & c& 0 \\ d & 0 & 0& a\end{pmatrix},
\label{R}
\ee
where $a, b, c, d$ are usually functions of spectral parameters $u$ and $v$.  The $\check{R}$-matrix can also be written as
\begin{align}
\check{R} & = \frac{1}{2} (a+c) \mathbbm{1} \otimes \mathbbm{1} +\frac{1}{2} (b+d)  \sigma^x  \otimes  \sigma^x \nonumber \\
&+\frac{1}{2} (b-d)  \sigma^y  \otimes  \sigma^y+ \frac{1}{2}( a- c) \sigma^z \otimes \sigma^z,
\label{R1}
\end{align}
which after a BST yields
\begin{align}
\check{G} & = \frac{1}{2} (a+ d \Sigma_2^x)( 1 + \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} (c+ b \Sigma_2^x)( 1- \Sigma_1^z \Sigma_3^z).
\label{G}
\end{align}
Similarly, the Kraus operators of Model B1 ($\tau \kappa =1$) are mapped to $K_0=  1/\sqrt{u+1} $ and
\begin{equation}
K_1=  \sqrt{\frac{u}{u+1}}[ \frac{1}{2} ( 1 + \Sigma_1^z \Sigma_3^z)+   \frac{1}{2}  \Sigma_2^x( 1- \Sigma_1^z \Sigma_3^z)].
\ee
Their action on the transformed $\rho$ is similar to Eq.(\ref{Kraus}). In the limit $u \to \infty$, the action is unitary with $K_1$ the same as RCA150 (see Appendix \ref{app1}). Other cases in Ref. \cite{de2021constructing} can be described similarly.

\section{Classical SCA}
\label{SCA}
In a standard cellular automaton\cite{Wolfram}, each site is one classical bit and it is updated according to a rule that takes the values of itself and its neighbors as the input \cite{farrelly2020review}. If the update rule is reversible, i.e., each output state has only one predecessor, it is also called a reversible cellular automaton; otherwise, it is irreversible.   Classical SCA can be obtained from QCA by constraining the density matrix $\rho$ to be diagonal, $\rho = \text{diag}(\rho_{\uparrow\uparrow}, \rho_{\uparrow\downarrow}, \rho_{\downarrow\uparrow}, \rho_{\downarrow\downarrow})$. That the restricted $\check{R}$ is stochastic (i.e., it maps a probability distribution to another one) is guaranteed by the CPTP property of the original $\check{R}$. (A unital $\check{R}$ reduces to a bi-stochastic matrix \cite{pavel2021correlations}.) Preservation of the diagonal space leads to classical (discrete) flows of probabilities $P(s_1,s_2,..., s_L) = \langle s_1, s_2,...,s_L| \rho |s_1, s_2,...,s_L\rangle$ and the continuum limit of the five circuits with two-site classical gates is described Ref.\cite{de2021constructing}. As an example, Model B1 with $\tau \kappa =1$, for finite $u$ yields
\begin{align}
\text{diag}&(\rho_{\uparrow\uparrow}, \rho_{\uparrow\downarrow}, \rho_{\downarrow\uparrow}, \rho_{\downarrow\downarrow}) \nonumber \\
&\to \text{diag}(\rho_{\uparrow\uparrow}, \frac{\rho_{\uparrow\downarrow} +u \rho_{\downarrow\uparrow}}{u+1}, \frac{\rho_{\downarrow\uparrow}+u \rho_{\uparrow\downarrow}}{u+1}, \rho_{\downarrow\downarrow}).
\end{align}
The corresponding transition matrix after a BST can be expressed as
\begin{align}
\check{G}(u) & = \frac{1}{2}( 1 + \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} \frac{1+ u \Sigma_2^x}{u+1}( 1- \Sigma_1^z \Sigma_3^z).
\end{align}

These integrable models are a subset of larger integrable family: ZFEV model (cf. Appendix \ref{app2}). Indeed, the application of a BST to the integrable chains Eq.(\ref{R1}) yields Eq.(\ref{G}).
The coefficients are given in Eq.(\ref{abcd}).
For $u \ge 0$ and $v \ge 0$, they are all nonnegative and $a + d = b + c =1$. They correspond to the ZFEV model with
\begin{align}
\Delta  &=  \frac{a^2 + b^2 - c^2 - d^2}{2(a b + c d)} = \text{sech} (2v), \nonumber \\
\Gamma & = \frac{a b - c d}{a b + c d} = \text{sech} (2v).
\end{align}
Since $\text{sech}(2v) \le 1$, the ground state of the ZFEV model is always in the disordered phase or on its critical boundaries \cite{baxter2016exactly}. In particular, RCA150 corresponds to $v=0$ after a BST, i.e., it always sits on a critical line. For a recent analysis, see Ref.\cite{wilkinson2021exact}. For the ZFEV model, there are many symmetries which preserve the integrability \cite{baxter2016exactly}. In particular, $a\leftrightarrow d$ and $b\leftrightarrow c$ maps RCA150 to RCA105. Note that in Ref.\cite{gombor2021integrable}, RCA150 and RCA105 were embedded into two different families of integrable quantum models. This is because the regularity condition was imposed. We show in Appendix \ref{app3} that the regularity condition can be relaxed without breaking the integrability and it corresponds exactly to the $a\leftrightarrow d$ and $b\leftrightarrow c$ symmetry. The transition matrix is mapped to
\begin{align}
\check{G}(u) & = \frac{1}{2}( 1 - \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} \frac{1+ u \Sigma_2^x}{u+1}( 1+ \Sigma_1^z \Sigma_3^z).
\end{align}

Thus, we can see that some integrable SCA can be embedded into a family of integrable nonunitary QCA. In particular, RCA150 (and RCA105) can be embedded into families of integrable SCA, unitary QCA \cite{gombor2021integrable}, and nonunitary QCA.

\section{Conclusion and outlook}
We have constructed several two-site translationally invariant integrable nonunitary quantum circuits using the $\check{R}$-matrices that appeared in constructing integrable Lindblad superoperators in the continuous case  in Ref. \cite{de2021constructing}. The trotterized Hubbard model with imaginary interaction strength \cite{sa2021integrable} is one special case of our construction. We showed that these $\check{R}$-matrices generate integrable open systems by demonstrating that they are CPTP, and thus represent a valid quantum channel. By applying the BST, integrable quantum circuits with three-site gates can also be constructed.  Among them, the classical RCA150 (and RCA105) can be embedded into an integrable family describing quantum open systems. The classical discrete flows on the diagonal of density matrices of the entire family of Model B1 can be embedded into the well-known integrable ZFEV model with parameter restricted such that it ground state is always in the disordered phase (and its boundary critical lines). Since the BST is a special case of the more general Clifford transformations \cite{jones2021integrable}, other descencant integrable nonunitary quantum circuits can be similarly derived.

In general, if we can find CPTP solutions to the YBE, we can use them directly to build integrable nonunitary quantum circuits. No Markovian approximation would be needed. Since we can regard the $\check{R}$-matrices of Model B1 and of Model B2 as those obtained from their unitary analogs via some analytical continuation, a closer look at the analytical continuation of known unitary solutions of YBE may may help find other CPTP solutions.

Our work points to a possible systematic way to obtain CPTP solutions to the YBE. Since all $\check{R}$-matrices of the five models in Ref.\cite{de2021constructing} are CPTP, a natural direction is to verify whether all $\check{R}$-matrices that can be constructed using the Lindbladian approach of Ref. \cite{de2021constructing} satisfy this property. In particular, if we use an ansatz of the form $\check{R}(u) =\mathbbm{1} + \sum_{n =1}^{\infty} r^{(n)}({\cal{L}}) u^n$ where $r^{(n)}({\cal{L}})$ is a polynomial of degree $n$ in ${\cal{L}}$ \cite{de2019classifying}, the trace-preserving property of $\check{R}$ is automatically guaranteed since $\text{Tr}[{\cal{L}}(\cdot)] =0$. Then it is sufficient to show that $\check{R}$ is completely positive.

A natural extension of the result in Ref.\cite{de2021constructing} to multiple families of jump operators may also lead to extension of the results in this paper to more general integrable open quantum circuits. In particular, it would be interesting to know if there is any connection to operator-space fragmentation discussed in Ref. \cite{Essler2020integrability} with two families of jump operators $l_{j, j+1}^{(1)} = \sqrt{J_1} \sigma_j^+\sigma_{j+1}^-$ and $l_{j, j+1}^{(2)} = \sqrt{J_2} \sigma_j^-\sigma_{j+1}^+$.

We have noticed that  both Model B1 and Model B2 have the Kraus representation $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$, a convex combination of unitary $U_i$. It superficially resembles the action of local gates in a random unitary quantum circuit where $U_i$  may be taken to be Haar-random unitary \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator}  or Clifford-random unitary \cite{li2018quantum, Ippoliti2021post}. Even though random systems are generically chaotic, it is tempting to see if techniques invented in the context of random circuits can be adopted to integrable nonunitary circuits. On the other hand, the action $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$ can be interpreted as an ensemble of quantum trajectories \cite{wiseman1996quantum, jian2020criticality}. If we focus on a single quantum trajectory and interpret nonunitary quantum gates as weak measurements following unitary evolution \cite{li2018quantum, Ippoliti2021post}, then we do not have to use vectorized density matrices. In particular, we can construct more integrable non-unitary quantum circuits by applying the spacetime duality to integrable unitary quantum circuits \cite{Ippoliti2021post}. All we need to is relax the inversion property (Eq.(\ref{inv})) and regularity condition (Eq.(\ref{regcon})), which can be easily done. We will leave it for future work.

We have only proved the integrability of the open quantum circuits constructed. How the extensive number of conserved charges dictate the thermalization in integrable nonunitary quantum circuits is an important question. Explicit computation of time correlation functions can be challenging since the gates are no longer unitary. It would be interesting to see if the Bethe ansatz approach can yield useful information about the complex spectrum of the circuits \cite{ziolkowska2020yang} and, moreover, how dominant complex eigenvalues controls the long-time dynamics \cite{claeys2022correlations}. We may draw some inspirations from solutions to integrable systems with boundary drives \cite{vanicat2018integrable}. Numerically study of spectral statistics as in Ref.\cite{sa2021integrable} may provide another perspective. It is also helpful to compute dynamical entanglement-related qualities, such as entanglement spread, entanglement negativity, of the models. It will certainly give us a better understanding of irreversible/nonunitary QCA \cite{farrelly2020review, brennen2003entanglement, piroli2020quantum}.

\section*{Acknowledgement}
We thank Aashish Clerk for discussions. This material is based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.

\section{Introduction}

Nonequilibruim quantum dynamics has been a focus of attention in theoretical and experimental physics in the recent years. Compared to their equilibrium counterparts, nonequilibrium systems are
more challenging to describe due to the explicit time-dependence of their evolution. That generally involves coherent superposition of a large number of energy eigenstates, and not only their statistical mixture.
One key question has to do with the ways an interacting quantum system reaches \cite{deutsch1991quantum, srednicki1994chaos} (or evades \cite{moudgalya2021quantum}) thermodynamic equilibrium.
Significant progress \cite{Polkovnikov2011,nandkishore2015many,d2016quantum, Abanin2019colloquium, Bertini2021finite, Alba2021genearlized} has been achieved in understanding of the thermalization both in continuous and discrete models, especially in one-dimensional (1D) quantum systems. Some advanced numerical methods have been developed \cite{schollwock2005density, schollwock2011density, Pekker2014hilbert, cirac2021matrix}, yet they remain limited in accessible system sizes and evolution times. A number of analytical results have been obtained, especially in integrable systems  where an extensive number of conserved quantities  strongly constrain the dynamics \cite{baxter2016exactly, vsamaj2013introduction, eckle2019models}. There, the generalized hydrodynamics approach, originally introduced to study inhomogeneous quantum quenches in integrable many-body systems, has been very fruitful \cite{Bertini2021finite, Alba2021genearlized}.   Exact analytical results have  been obtained recently in systems including the Lieb-Liniger model \cite{Granet2021systematic} and the folded XXZ model \cite{pozsgay2021integrable, Zadnik2021folded1, Zadnik2021folded2}.

Typically, integrating time evolution involves discretization of time, with the small time step  ensuring that only small changes in the state of the system during a time step occur.
However, discretizing time, without assuming smallness of the time step, can also be useful for modeling physical systems.
Thus obtained quantum cellular automata (QCA) are locality-preserving maps that act in discrete space and time \cite{farrelly2020review, piroli2020quantum}, and can be regarded as a special kind of tensor networks \cite{piroli2020quantum, cirac2021matrix}. The focus of this paper are QCA that are constant-depth quantum circuits, which consist of layers of gates acting on a lattice of qubits \cite{farrelly2020review, cirac2021matrix}. Several exact results have already been obtained in such quantum circuits, including random quantum circuits \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator}, dual unitary circuits \cite{bertini2019exact, bertini2019entanglement, piroli2020exact, pavel2021correlations, claeys2021ergodic} and the Rule 54 quantum cellular automaton (QCA54)\cite{friedman2019integrable, klobas2021exact, klobas2021exact2, klobas2021ent}. The first two types are generically chaotic while QCA54 is believed to be integrable \footnote{However, it is not clear if the model can fit into a family of Yang-Baxter integrable QCA (see Ref. \cite{gombor2021integrable, prosen2021reversible} for more details).}. In QCA54, relaxation dynamics and entanglement dynamics after a quantum quench from some initial states can be computed explicitly \cite{klobas2021exact, klobas2021exact2, klobas2021ent}.

Integrable circuits can be constructed by using $\check{R}$-matrices that satisfy the Yang-Baxter equation (YBE) as the local gates  (Fig.~\ref{fig1}).  Depending on the properties of the $\check{R}$-matrix, there are three possible cases: (i) classical stochastic cellular automata (SCA) if $\check{R}$ is stochastic; (ii) unitary QCA if $\check{R}$ is unitary; (iii) nonunitary QCA if
$\check{R}$ is completely positive and trace preserving (CPTP). Most of the previous work focused on the first and the second case, while integrable nonunitary QCA are less studied.

Nonunitary quantum circuits \cite{farrelly2020review, brennen2003entanglement, piroli2020quantum}, especially unitary-projective ones \cite{li2018quantum, skinner2019meas, chan2019unitary, jian2020criticality, Ippoliti2021post, lu2021spacetime, Ippoliti2022fractal}, are important not only because they correspond to discrete-time versions of open quantum systems, but also since they can represent gate-based quantum computing and  can be used to study error-mitigation algorithms in near-term quantum devices.
One example of an integrable nonunitary QCA was recently obtained by trotterizing the Hubbard model with imaginary interaction strength \cite{sa2021integrable}. Such examples are highly nontrivial to build, which leaves a question whether there is a systematic way to construct integrable nonunitary quantum circuits.

In this paper, we construct several more examples of integrable nonunitary quantum circuits by making use of the recent developments in the systematic construction of Yang-Baxter(YB) integrable Lindblad superoperators in the continuous-time case \cite{ziolkowska2020yang, de2021constructing}. The procedure is based on finding $\check{R}$-matrices that satisfy the YBE, such that their derivatives at special points correspond to integrable Lindblad superoperators. What we demonstrate here is that thus obtained $\check{R}$-matrices not only yield continuous modes, but can also be used as the building blocks of discrete-time nonunitary quantum circtuits.
In combination with the approach of Ref.  \cite{de2021constructing} this may pave the way for a more systematic construction of nonunitary integrable quantum circuits.

The rest of the paper is organized as follows: In Section \ref{sect2}, we review the construction of integrable nonunitary quantum circuits with two-site gates and continuous integrable Lindblad superoperators. By demonstrating the CPTP property of the $\check{R}$-matrices constructed from the latter, we show that these  $\check{R}$-matrices can be used to build integrable nonunitary quantum circuits. In Section \ref{sect3}, we apply the bond-site transformation to the quantum circuits with two-site gates to construct integrable nonunitary quantum circuits with three-site gates. In Section \ref{SCA}, we briefly discuss classical SCA when the density matrix is restricted to be diagonal. In particular, we show that the Rule 150 reversible cellular automaton (RCA150) \cite{bobenko1993two} can be embedded into one of the models discussed and the Rule 105 reversible cellular automaton (RCA105) can also be embedded if the regularity condition is relaxed. In the last section, we summarize and discuss a potential way to build integrable nonunitary quantum circuits systematically, and point the way for possible future directions. Some details are relegated to the appendixes.

\section{Integrable quantum circuits with two-site gates}
\label{sect2}
\subsection{Integrability structure}
\label{int_str}
In this section, we review the integrability structure of (two-site) translationally invariant quantum circuits with two-site gates \cite{vanicat2018integrable, sa2021integrable}. It is based on the standard quantum inverse-scattering method
\cite{vsamaj2013introduction, eckle2019models, baxter2016exactly}. The central idea is to find a family of mutually commuting transfer matrices obtained from an $\check{R}$-matrix that satisfies the braid YBE. Then we can obtain a family of commuting conserved charges and even use algebraic Bethe ansatz techniques to obtain the spectrum of the transfer matrices and eigenfunctions.  For a quick introduction, see, e.g., Ref. \cite{Bertini2021finite}.

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{circuit2}
\caption{Floquet quantum circuit with period two. Each site has a Hilbert space of dimension $d$. Two-site gates $\check{R}$ that satisfy the YBE (Eq.(\ref{YBE})) act alternately.}
\label{fig1}
\end{figure}

Assume that we have an $\check{R}$-matrix, parametrized by $u$ and $v$, acting on two adjacent sites with local Hilbert space dimension $d$ each and that the $\check{R}$-matrix satisfies the braid YBE \cite{vsamaj2013introduction}:

\begin{align}
\check{R}_{23}(u, v)& \check{R}_{12}(u,w)\check{R}_{23}(v,w) =\nonumber \\
&\check{R}_{12}(v, w)\check{R}_{23}(u,w)\check{R}_{12}(u,v).
\label{YBE}
\end{align}
We let it act on each pair of sites $2j-1$ and $2j$ denoted by $(2j-1, 2j)$ with $j \in\{1, 2,..., L/2\}$ of a lattice of even length $L$ (Fig.~\ref{fig1}). We take periodic boundary conditions, namely, we identify $2L+1\sim 1$. Then we let it act on all pairs $(2j, 2j+1)$ of sites for $j \in\{1, 2,..., L/2\}$. The combination of the actions is given by
\begin{equation}
{\cal{V}}(u_1, u_2)  = \prod_{j=1}^{L/2}\check{R}_{2j, 2j+1} (u_1, u_2) \prod_{j=1}^{L/2}\check{R}_{2j-1, 2j} (u_1, u_2).
\ee
Then we repeat the action. In this way, we construct a brickwork type Floquet  quantum circuit with period 2. If we take $d =2$ for a closed system, then $\check{R}$ is 4 by 4. If we are to study open systems, it is natural to use density matrix formalism. Namely, the Floquet operator ${\cal{V}}$ acts on density matrix $\rho$. After vectorizing $\rho$, $\rho = \sum \rho_{ij} |i\rangle \langle j| \to |\rho\rangle=   \sum \rho_{ij} |i\rangle | j\rangle\rangle$, the local dimension becomes $d \to d^2$, and thus $\check{R}$ is 16 by 16.  In this section, we impose the regularity condition
\begin{equation} \check{R}(u, u) = \mathbbm{1}, \label{regcon}\ee
or equivalently, \begin{equation} R(u, u) \equiv P \check{R}(u, u) =P,\ee where $P$ is the permutation operator between two spaces defined by $P|i\rangle |j\rangle = |j\rangle |i\rangle$, which is also called the SWAP gate. It describes the trivial exchange process of the scattering of two particles \cite{vsamaj2013introduction}.  This condition, together with the YBE, implies \cite{gombor2021integrable} the inversion property,
\begin{equation} \check{R}_{12}(u, v) \check{R}_{12}(v, u) =\mathbbm{1}.\label{inv}\ee

To show that the circuit is YB integrable, we follow Ref.~\cite{sa2021integrable} and define a staggered transfer matrix:
\begin{align}
T(u) & =  \nonumber \\
&\text{Tr}_{a}[R_{a L} (u, u_1 )R_{a L-1} (u, u_2 )... R_{a 2} (u, u_1 ) R_{a 1} (u, u_2 )]
\end{align}
where $a$ can be interpreted as an auxiliary space to take trace over. Then the YBE leads to
$[T(u), T(v)]=0$ for any $u$ and $v$.

The regularity property and the inversion property then imply
\begin{align}
 &T(u_1)   \nonumber \\
=&P_{12} P_{13}...P_{1L}   \check{R}_{12} (u_1, u_2 )\check{R}_{34} (u_1, u_2 )... \check{R}_{L-1 L} (u_1, u_2 ) \nonumber \\
 \equiv & U \check{R}_{12} (u_1, u_2 )\check{R}_{34} (u_1, u_2 )... \check{R}_{L-1 L} (u_1, u_2 ),
\end{align}
and
\begin{align}
& T(u_2) \nonumber \\
=& P_{12} P_{13}...P_{1L} \check{R}_{23}^{-1} (u_1, u_2 )\check{R}_{45}^{-1} (u_1, u_2 )... \check{R}^{-1}_{L 1} (u_1, u_2 ) \nonumber \\
 \equiv & U  \check{R}_{23}^{-1} (u_1, u_2 )\check{R}_{45}^{-1} (u_1, u_2 )... \check{R}^{-1}_{L 1} (u_1, u_2 ),
\end{align}
where $U =P_{12} P_{13}...P_{1L}$ is a one-site shift. We have used the identity ${\cal{O}}_{ab} P_{ac} = P_{ac} {\cal{O}}_{cb}$ (for $b\neq c$) and $\text{Tr}_a P_{a1} = \mathbbm{1}$. Thus,

\begin{equation}
{\cal{V}} (u_1, u_2)  = T(u_2)^{-1} T(u_1)
\ee
and it commutes with $T(u)$ for any $u$. The commuting charges are given by
\begin{equation}
Q_n^+ =\frac{d^n}{d u^n}  \ln T(u) |_{u =u_1}, \quad Q_n^- =\frac{d^n}{d u^n}  \ln T(u) |_{u =u_2}.
\ee

\subsection{Integrable Lindblad Superoperators}\label{ILS}
Quantum circuits that we discussed may be regarded as coming from trotterization of  corresponding lattice models. If an open continuum system is in contact with a Markovian environment whose dynamics is much faster than  the internal one, the equation of motion is approximately given by the Lindbladian master equation
\begin{equation}
\dot{\rho} = {\cal{L}} \rho = i[\rho, H] + \sum_a l_a \rho l_a^{\dagger} - \frac{1}{2}\{ l_a^{\dagger} l_a, \rho\},
\label{lme}
\ee
where $\rho$ is the density operator, ${\cal{L}}$ is the Lindbladian and $l_a$ are jump operators. For a lattice system,  ${\cal{L}} = \sum_i {\cal{L}}_{i, i+1}$.

A recent study \cite{de2021constructing} constructed many new examples of integrable Lindblad superoperators following Ref.~\cite{de2019classifying, de2020classifying}.  We briefly summarize their approach. In terms of the vectorized density matrix, $\rho \in {\cal{H}}^{(1)} \otimes  {\cal{H}}^{(2)} $,  the Lindbladian density for a system with one family of jump operators has the explicit form
\begin{align}
{\cal{L}}_{i, i+1}(u) & = -i h_{i, i+1}^{(1)}(u) + i h_{i, i+1}^{(2)*}(u) +l_{i, i+1}^{(1)}(u) l_{i, i+1}^{(2)*}(u)  \nonumber \\
& - \frac{1}{2} l_{i, i+1}^{(1)\dagger}(u) l_{i, i+1}^{(1)}(u) - \frac{1}{2} l_{i, i+1}^{(2)\text{T}}(u) l_{i, i+1}^{(2)*}(u).
\label{lin}
\end{align}
Here, ${\cal{H}}^{(1)}$ is that of two qubits and ${\cal{H}}^{(2)}$ is its canonical dual space; $h^{(s)}$ and $l^{(s)}$ act on the $s$-th Hilbert space. $^*$ denotes the complex conjugate and $^\text{T}$ transpose. The dependence on the spectral parameter $u$ will be used in the following.

To construct integrable ${\cal{L}}$, we use ansatze for ${\cal{L}}$ of the form in Eq.(\ref{lin}) and its associated $\check{R}$ satisfying
\begin{equation}
\check{R}(u, u) =\mathbbm{1}, \quad   {\cal{L}}_{i, i+1}(u) = \partial_v \check{R}_{i, i+1}(v, u)|_{v\to u}.
\label{initial}
\ee

Let $Q_n$ be the $n$-site conserved charge, in particular, $Q_2 = {\cal{L}}$. Higher charges are generated by
\begin{equation}
Q_{n+1}(u) = [ B(u), Q_n(u)], \quad n >1,
\label{q2}
\ee
where the boost operator $B(u):=  \partial_u + \sum_{i = -\infty}^{\infty} i {\cal{L}}_{i, i+1}(u)$. This expression is defined for infinite chains but reduces to finite chains consistently. Then impose the condition
\begin{equation} [Q_2(u), Q_3(u)] =0\ee on the ansatz and solve the following equations obtained from the YBE:
\begin{subequations}
\begin{gather}
[R_{13} R_{23}, {\cal{L}}_{12}(u)] = \dot{R}_{13} R_{23} - R_{13} \dot{R}_{23};\\
[R_{13} R_{12}, {\cal{L}}_{23}(v)] = R_{13} R'_{12} - R'_{13} R_{12}.
\end{gather}
\label{suther}
\end{subequations}
In the first equation $R_{ij} = R_{ij}(u, u)$ and $\dot{R}$ is the derivative taken with respect to the first varaible at $v=u$. In the second equation, $R_{ij} = R_{ij}(v, v)$ and the derivative is taken with respect to the second variable at $u =v$. Combining Eqs.(\ref{initial}), (\ref{q2}) and (\ref{suther}), we can obtain $\cal{L}$ and $\check{R}$. In the end, we can check if the $\check{R}$-matrix obtained this way satisfies the YBE (Eq.(\ref{YBE})).

Note that there are some degrees of freedom in $h$ and $l$ that keep ${\cal{L}}_{i, i+1}(u)$ invariant \cite{de2021constructing}, thus yielding the same $\check{R}$. However, different ${\cal{L}}_{i, i+1}(u)$ may give rise to $\check{R}$-matrices that are still equivalent via local basis transformations, reparameterizations, normalizations, discrete transformations, and twists. See Ref.\cite{de2019classifying, de2020classifying} for more information. In order to qualify as integrable non-unitary (open) quantum circuits,  the $\check{R}$-matrix has to be a quantum channel or, equivalently, satisfy CPTP \cite{preskill1998lecture}.

\subsection{CPTP}
To show that an $\check{R}$-matrix is CPTP, we can use the Choi-Jamiolkowski representation \cite{sa2021integrable} of $\check{R}$ and define the Choi matrix $D_{\beta \delta \zeta \theta}^{\alpha \gamma\varepsilon \eta} \equiv \check{R}_{\gamma \delta \eta \theta}^{\alpha \beta \varepsilon \zeta}$. Then, $\check{R}$ is trace-preserving if and only if  $D^{\alpha \gamma \varepsilon\eta}_{\alpha \delta \varepsilon \theta} = \check{R}^{\alpha\alpha \varepsilon\varepsilon}_{\gamma \delta \eta \theta} = \delta_{\alpha}^{\gamma} \delta_{\theta}^{\eta}$; $\check{R}$ is completely positive if and only if $D$ is non-negative. Additionally,  $\check{R}$ is unital if and only if $D^{\alpha \gamma \varepsilon\eta}_{\beta \gamma \zeta \eta} = \check{R}^{\alpha\beta \varepsilon\zeta }_{\gamma \gamma \eta \eta} = \delta^{\alpha}_{\beta} \delta^{\varepsilon}_{\zeta}$. Unitality is not necessary for $\check{R}$ to be CPTP.

It is easy to check explicitly that all $\check{R}$-matrices of all five cases in Ref.~\cite{de2021constructing} are CPTP. Additionally,  \begin{equation} \text{tr}(D (u))= \text{tr}(D(0)) =4.
\ee
Hence, it has at most rank 2 and the two non-vanishing eigenvalues satisfy
\begin{equation}
\lambda_1 +\lambda_2 =4,
\ee
and $\lambda_1, \lambda_2 \ge 0$ for $u\ge 0$. Since $\check{R}$ is CPTP, we have a Kraus representation \cite{preskill1998lecture} of  $\phi_{\beta \delta \zeta \theta}^{\alpha \gamma\varepsilon \eta} \equiv \check{R}_{\beta \zeta \delta \theta}^{\alpha \varepsilon \gamma \eta}$:
\begin{equation}
\phi  =K_0\otimes K_0^* + K_1\otimes K_1^* ,
\ee
which satisfies
\begin{equation}
K_0^{\dagger}K_0 +K_1^{\dagger}K_1 = \mathbbm{1}.
\ee
It acts on $\rho$ as
\begin{equation}
\phi[\rho] = K_0\rho K_0^{\dagger} + K_1 \rho K_1^{\dagger}.
\ee
$\check{R}$ is unital if
\begin{equation}
K_0 K_0^{\dagger} +K_1  K_1^{\dagger}= \mathbbm{1}.
\label{unital}
\ee

Thus, we can use these $\check{R}$-matrices to construct integrable nonunitary quantum circuits.
Note that the CPTP property of $\check{R}$ holds for all $u \ge 0$. However, when $u\to 0$, the quantum circuits reduce to the continuous systems described by the Lindbladian master equation (Eq.(\ref{lme})). Indeed, in a closed system with $d=2$, the $u\to 0$ limit of the quantum circuit described in Subsection \ref{int_str} reduces to $\lim_{N\to \infty} {\cal{V}}(u/N)^N  \to \exp(i u H)$, where $H$ can be interpreted as the hermitian Hamiltonian \cite{vanicat2018integrable}. Similarly, in the open system, when we take $u\to 0$, $ \lim_{N\to \infty} {\cal{V}}(u/N)^N \to \exp(u{\cal{L}})$.
The quantum circuit reduces to a system with continuous `time' $u$ \cite{vanicat2018integrable}.

In general, if there are more than one family of jump operators, the expansion of the Kraus operators around $u =0$ is  \footnote{Note that if $\check{R}$ has a difference form, i.e., $\check{R}(u, v) = \check{R}(u-v)$, we simply write $\check{R}(u) = \check{R}(u, 0)$. Similarly for $K_a$ and other operators. In this case, ${\cal{L}}$ is independent on $u$.}
\begin{subequations}
    \begin{align}
K_0(u) & = I + (-i H - \frac{1}{2} \sum_{a >0} l_a^{\dagger} l_a ) u + O(u^2), \\
K_a (u) & = \sqrt{u} l_a  + O(u).
\end{align}
\label{kraus}
\end{subequations}

In the following we show that the five $\check{R}(u)$-matrices in  Ref.\cite{de2021constructing} constructed using integrable Lindblad superoperators are CPTP. Indeed, it would be interesting to see if all $\check{R}$-matrices constructed in this way are always CPTP; so far we have not been able to demonstrate that.
In principle, we can also run the argument in reverse: If we have a YBE solution $\check{R}(u)$ that has the CPTP property, it can be used to obtain a corresponding integrable Lindblad superoperator.

\subsection{Examples}
We list all five cases in
Ref.~\cite{de2021constructing}. We can classify them into two categories. The first category includes Model B1 and Model B2 where two limiting cases ($u\to 0$ and $u\to \infty$) are both unitary. The rest (Model A1, Model A2, and Model B3) do not have this property.
\subsubsection{Model B1}
The (16 by 16) $R$-matrix we use is obtained from that given in Ref.~\cite{de2021constructing} by rescaling it by a factor $e^u/(u+1)$. $h =0$ and
\begin{equation}
l =
\begin{pmatrix} \tau & 0 & 0 & 0\\
0 & 0& 1 & 0 \\
0 &1 &0  & 0\\
0 & 0 & 0 & \kappa
\end{pmatrix},
\ee
where $\tau = \pm 1$ and $\kappa =\pm 1$. For $\tau =\kappa =1$,
\begin{equation}
\check{R}(u) = \frac{\mathbbm{1} +u P}{u+1},
\ee
where $P$ is the permutation matrix.
This is the $R$-matrix of the $SU(4)$ invariant chain if $u$ is taken to be purely imaginary. It is easy to see that $\check{R} = P R$ satisfies the inversion property Eq.(\ref{inv}).

In this case, the $\check{R}$-matrix is both CPTP and unital. The non-trivial eigenvalues of the Choi matrix $D$ are
\begin{equation}
\lambda_{\pm} = 2 \pm 2 \frac{\sqrt{u^2 -(\frac{3-\tau \kappa}{2} )u+1}}{u+1}.
\ee

The Kraus operators are
\begin{equation} K_0(u) = \sqrt{\frac{1}{u+1}} \mathbbm{1}_{4\times 4}, K_1(u) = \sqrt{\frac{u}{u+1}}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& \tau & 0 \\
0 &\tau &0  & 0\\
0 & 0 & 0 & \kappa\tau
\end{pmatrix}.
\ee

The action on $\rho$ is
\begin{equation}
\phi[\rho] = \frac{\rho + u K_1 \rho K_1 }{u+1}.
\label{Kraus}
\ee
It is easy to check that  $K^{\dagger}_0 K_0 + K_1^{\dagger} K_1 =\mathbbm{1}$.  When $u =0$, $K_0(0) = \mathbbm{1}$ and $K_1(0) =0$. When $u\to \infty$,
\begin{equation} K_0(\infty) = 0, \quad K_1(\infty) =
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& \tau & 0 \\
0 &\tau &0  & 0\\
0 & 0 & 0 & \kappa\tau
\end{pmatrix}.
\label{k0k1}
\ee
Thus, the trivial unitary quantum circuit $\check{R} =\mathbbm{1}$ is connected to the unitary circuit $U = K_1(\infty)$ by tuning $u$ from 0 to $\infty$.  For $\tau \kappa =1$, $K_1(\infty)$ corresponds to the XXZ model with $\Delta =\pm 1$. In this case, $K_1(\infty)$ will not generate entanglement for any (pure) product states, which is not true for the other two cases. For $\tau \kappa =- 1$, it corresponds to the XX model with a specific $h_z$. For all four cases, $K_1(\infty)$ is dual-unitary, and it belongs to the noninteracting case described in Ref. \cite{bertini2019exact}.

That $\check{R}$ is unital implies that the maximally mixed state is a steady state. It is the quantum version of the bi-stochastic classical map \cite{pavel2021correlations}. The conserved charges can be obtained as in Ref.\cite{vanicat2018integrable}.

\subsubsection{Model B2}
This case includes the XX model with dephasing noise (which corresponds to the Hubbard model with imaginary coupling) discussed in Ref.~\cite{sa2021integrable} as the special case ($\phi =\pi/2$ in the following discussion). $h$ and $l$, respectively, are given below:
\begin{subequations}
\begin{equation}
h =2  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\begin{align}
&\frac{l}{\eta} = \nonumber \\
& \begin{pmatrix} \cosh(2u) & 0 & 0 & 0\\
0 & 1& \text{i} e^{\text{i} \phi} \sinh(2u) & 0 \\
0 & -\text{i} e^{-\text{i} \phi} \sinh(2u)& -1  & 0\\
0 & 0 & 0 & -\cosh(2u)
\end{pmatrix},
\end{align}
\end{subequations}
where \begin{equation} \eta (u) = \sqrt[4]{\frac{w^2}{ w^2\sinh^2 (2u) +1}}.\ee

Define
\begin{equation}
\check{G}(u) =\frac{1}{a}  \begin{pmatrix} a & 0 & 0& 0\\ 0 & c &-\text{i} b e^{\text{i} \phi} & 0 \\  0 & -\text{i} b e^{-\text{i} \phi} & c& 0 \\ 0 & 0 & 0& a\end{pmatrix},
\ee
where \begin{equation} a = \cosh (u), \quad b = \sinh(u), \quad c=1.\ee
Also, take the matrix $e^{\beta}_{\alpha}$ whose matrix element is $(e^{\beta}_{\alpha})^{\beta'}_{\alpha'} =\delta_{\beta, \beta'}\delta_{\alpha, \alpha'}$, and
\begin{subequations}
\begin{gather}
\check{r}_{\uparrow}(\lambda) = \check{G}^{\alpha \gamma}_{\beta \delta}(\lambda) e_{\alpha}^{\beta} \otimes \mathbbm{1}_2 \otimes e_{\gamma}^{\delta} \otimes \mathbbm{1}_2, \\
\check{r}_{\downarrow}(\lambda) = \check{G}^{\alpha \gamma}_{\beta \delta}(\lambda)\mathbbm{1}_2\otimes e_{\alpha}^{\beta} \otimes \mathbbm{1}_2 \otimes e_{\gamma}^{\delta},
\end{gather}
\end{subequations}
Define $\check{r} = \check{r}_{\uparrow} \check{r}^*_{\downarrow}$ and then the $\check{R}$-matrix can take the form
\begin{equation}
\check{R}(u, v) = \beta \check{r}(u-v) + \alpha \check{r}(u+v) (\sigma^z \otimes \sigma^z \otimes \mathbbm{1}_4),
\ee
where
\begin{equation}
\beta =
  \frac{\cosh(u - v) \cosh(
     \xi - \eta)}{(\cosh (u - v) \cosh( \xi - \eta) + \cosh(u + v) \sinh( \xi - \eta))},
\ee
$\alpha = 1-\beta$, and
\begin{equation}
\frac{\sinh(2\xi)}{\sinh(2u)} =\frac{\sinh(2\eta)}{\sinh(2v)} = w, \quad u\ge v.
\ee
The inversion condition, Eq.(\ref{inv}), is satisfied. $\check{R}(u, v)$ is not of difference form  and it is unital and CPTP. In particular, the nonvanishing eigenvalues of its Choi matrix are $4\alpha$ and $4\beta$. The Kraus operators are given by
\begin{equation}
K_0 = \sqrt{\beta}   \check{G}(u-v),\quad  K_1 = \sqrt{\alpha}   \check{G}(u+v) (\sigma^z \otimes \mathbbm{1}).
\ee
When $u=v$, $\xi =\eta$, $\beta =1$ and $\alpha=0$. $K_0(0)=\mathbbm{1}$ and $K_1(0) =0$. When $u \to \infty$, $\alpha =1$ and $\beta=0$. Then
\begin{equation} K_0(\infty) = 0, \quad K_1(\infty) =
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& -i e^{i \phi} & 0 \\
0 &i e^{-i \phi} &0  & 0\\
0 & 0 & 0 & -1
\end{pmatrix}.
\ee
Note that when $\phi =\pi/2$ this is same as Eq.(\ref{k0k1}) for $\tau=-1$ and $\kappa =-1$. Both limiting cases are unitary. The integrability of the circuit is also demonstrated numerically in Ref.\cite{sa2021integrable} by examining the spectral statistics.

Note that both Model B1 and Model B2 have the following form
$\phi[\rho] = p U_0\rho U_10^{\dagger} + (1-p) U_1\rho U_1^{\dagger}$ with $0\le p \le 1$.  More generally, $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$ has the form of a convex combination where $U_i$ is unitary. It resembles random unitary quantum circuits if we take $U_i$ to be a Haar-random unitary \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator} or  Clifford-random unitary \cite{li2018quantum, Ippoliti2021post}.

\subsubsection{Model A1}
In this model, $h$ and $l$ are given by
\begin{equation}
h =\frac{1}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},\quad
l=\\
 \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& 0& 0 \\
0 & 1& -i e^{-\text{i} \phi}  & 0\\
0 & 0 & 0 &0
\end{pmatrix}.
\ee
The entries of the $R$-matrix are given in Ref. \cite{de2021constructing}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital in this case since Eq.(\ref{unital}) is not satisfied. The two nonvanishing eigenvalues of the Choi matrix are $
\lambda_{\pm} =2 \pm \sqrt{1 +2 e^{-2 u} + e^{-3 u}}$. Thus
$\lambda_{\pm} \ge 0$ when $u \ge 0$.
The Kraus operators are given by
\begin{subequations}
    \begin{align}
 K_0 = e^{-u/2}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & e^{-u/2} & 0 & 0 \\
0 &0 &e^{u/2}  & 0\\
0 & 0 & 0 & 1
\end{pmatrix}, \\
 K_1 =  \sqrt{1-e^{-u}}
\begin{pmatrix} 1& 0 & 0 & 0\\
0 & e^{-u/2} & 0 & 0 \\
0 &-\text{i} e^{ \text{i} \phi} & 0 & 0\\
0 & 0 & 0 & 1
\end{pmatrix}.
\end{align}
\end{subequations}
The phases of $K_0$ and $K_1$ are irrelevant.

\subsubsection{Model A2}
Here, $h$ and $l$ are given by
\begin{equation}
h =\frac{\tau}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& i  & 0 \\
0 &-i &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},\quad
l=\\
 \begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0& 0& 0 \\
0 & \tau & 1  & 0\\
1 & 0 & 0 &0
\end{pmatrix},
\ee
where $\tau =\pm 1$.
The entries of the $R$-matrix of the Choi matrix are given in Ref. \cite{de2021constructing}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital in this case since Eq.(\ref{unital}) is not satisfied.
The two nonvanishing eigenvalues are
$\lambda_{\pm} = 2 \pm 2 e^{-3u/2}$ and $\lambda_{\pm} \ge 0$ for $u \ge 0$. The Kraus operators are given by
\begin{subequations}
    \begin{align}
    K_0= e^{-u}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & e^{u/2} & 0 & 0 \\
0 &-\tau (e^u-1)  &e^{u/2}  & 0\\
0 & 0 & 0 & e^{u}
\end{pmatrix}, \\ K_1 = \sqrt{e^u -1} e^{-u}
\begin{pmatrix} 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
0 &\tau &e^{u/2}  & 0\\
e^{u/2} & 0 & 0 & 0
\end{pmatrix}.
\end{align}
\end{subequations}
 $\check{R}$ is not unital in this case since Eq.(\ref{unital}) is not satisfied.

\subsubsection{Model B3}
$h$ and $l$ are given by
\begin{subequations}
 \begin{equation}
h  =\frac{1}{2}  \begin{pmatrix} 0 & 0 & 0 & 0\\
0 & 0& e^{\text{i} \phi} & 0 \\
0 &e^{-\text{i} \phi} &0  & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\begin{equation}
l =\sqrt{\frac{\gamma}{2}}
 \begin{pmatrix} \gamma & 0 & 0 & 0\\
0 & 1& i (\gamma-1) e^{i \phi}& 0 \\
0 & -i (\gamma+1) e^{-i \phi}& -1  & 0\\
0 & 0 & 0 &\gamma
\end{pmatrix}.
\ee
\end{subequations}
$\gamma\ge 0$.
The entries of the $R$-matrix of the Choi matrix are given in Ref. \cite{de2021constructing}   \footnote{There are typos for the $R$-matrix of Model B3 in Ref. \cite{de2021constructing} : $R_4^{10} = e^{\text{i} \phi} 2\text{i} (\gamma+1)^2 \sinh(\alpha)/\zeta$ and $R_4^7 = -e^{-\text{i} \phi} 2\text{i} (\gamma+1)^2 \sinh(\alpha)/\zeta$.}. $\check{R}$ is CPTP  and satisfies the inversion condition, but is not unital for generic $\gamma$ in this case since Eq.(\ref{unital}) is not satisfied.  The Kraus operators are given by

\begin{subequations}
    \begin{align}
   &  K_0 = \nonumber \\
     & \frac{1}{\eta}
\begin{pmatrix} \eta  & 0 & 0 & 0\\
0 & 1 & -\text{i} (1-\gamma)\sinh (\alpha)e^{i \phi} & 0 \\
0 & -\text{i} (1+\gamma)\sinh (\alpha) e^{-i \phi} & 1 & 0\\
0 & 0 & 0 &  \eta
\end{pmatrix},
\end{align}
\begin{equation}  K_1 = \frac{\xi}{\eta}
\begin{pmatrix} 0 & 0 & 0 & 0\\
0 & -e^{-\alpha} (1-\gamma) & i(1-\gamma) e^{i \phi}& 0 \\
0 & i(1+\gamma)e^{-i \phi} &  e^{\alpha} (1+\gamma) & 0\\
0 & 0 & 0 & 0
\end{pmatrix},
\ee
\end{subequations}
where $\alpha  = (\gamma^2+1)u/2$, $\eta  = \gamma \sinh(\alpha) +\cosh(\alpha)$, and $\xi= \sqrt{\gamma \sinh(\alpha) /[(1+\gamma^2) \cosh(\alpha) +2\gamma \sinh(\alpha)] }$.

\section{Integrable quantum circuits with three-site gates}
\label{sect3}

In the previous section we showed how two-site $\check{R}$ matrices that appear in the integrable Lindbladian construction can be repurposed to build non nonunitary quantum circuits. In this section, we extend this construction to integrable three-site nonunitary gates and circtuits by taking advantage of bond-site transformations.

\subsection{Bond-site transformation}
\label{sec_bst}
A bond-site transformation (BST) maps a relative orientation of a pair of spin to a single  spin value.  We discuss BST following Ref.~\cite{gombor2021integrable}. For a 1D spin chain, we choose the computational basis $\mid s_1, s_2, ..., s_L\rangle$ where $s_i =\pm 1$. For each pair of neighboring spins, define a state $s'_i = s_i s_{i+1}$ on the bond. More concretely, if two neighboring sites are in the same state, then $s' =1$; if two neighbor sites are in different states,then  $s' =-1$. Thus, the state $s'=-1$ can be interpreted as the existence of a domain wall. Note that this is a two-to-one correspondence: states related by a global spin reflection are mapped to the same state on the bond. If open boundary conditions are imposed on the original spin chain and one of the boundary states is fixed, the states on the bonds after the BST are in the one-to-one correspondence to the rest of the states. If periodic boundary conditions are taken, the mapping is two-to-one, and the bond states are constrained by $s'_1 s'_2\cdots s'_L =1$. Note that, similar to the Jordan-Wigner transformation \cite{fradkin2013field}, a BST is nonlocal because a domain wall is nonlocal in the original basis. Moreover, an entangled state can be mapped to an unentangled one. For example, the Greenberger–Horne–Zeilinger state $(\mid\uparrow\rangle^{\otimes L} +\mid\downarrow\rangle^{\otimes L})/\sqrt{2}$ is mapped to a product state.

On the operator level, an $(l+1)$-site operator may be mapped to an $l$-site operator after a BST. In particular, three-site operators that commute with the global spin reflection ($\Sigma^x_1\otimes \Sigma^x_2 \otimes \Sigma^x_3$) are in one-to-one correspondence with a two-site operators that commute with $\sigma^z \otimes \sigma^z$.

\begin{figure}[H]
\centering
\includegraphics[width=0.2\textwidth]{b.s.t.1}
\label{figbts}
\end{figure}

They are listed below:
\begin{align}
\sigma^x & \otimes \sigma^x  \leftrightarrow \mathbbm{1} \ \otimes \Sigma^x \otimes\ \mathbbm{1}, \quad
\sigma^y \otimes \sigma^y  \leftrightarrow - \Sigma^z \otimes \Sigma^x \otimes \Sigma^z,\nonumber  \\
\sigma^z& \otimes \sigma^z  \leftrightarrow \Sigma^z \otimes\ \mathbbm{1} \ \otimes \Sigma^z, \quad
\sigma^y   \otimes \sigma^x  \leftrightarrow \Sigma^z \otimes \Sigma^y \otimes\ \mathbbm{1},\nonumber \\
\sigma^x &\otimes \sigma^y  \leftrightarrow\ \mathbbm{1} \  \otimes \Sigma^y \otimes \Sigma^z, \quad
\mathbbm{1}   \otimes \sigma^z  \leftrightarrow -\mathbbm{1}\ \otimes \Sigma^z \otimes \Sigma^z,\nonumber \\
\sigma^z &\otimes \mathbbm{1} \leftrightarrow -\Sigma^z \otimes \Sigma^z \otimes\ \mathbbm{1}.
\end{align}
Consequently, general operators of the following form
\begin{align}
\check{R} &= h_1 \mathbbm{1} \otimes \mathbbm{1} + \frac{h_2}{2} (\sigma^z \otimes \mathbbm{1} - \mathbbm{1} \otimes \sigma^z) + h_3 \sigma^+ \otimes \sigma^- \nonumber \\
&+ h_4 \sigma^- \otimes \sigma^+ + \frac{h_5}{2} (\sigma^z \otimes \mathbbm{1} + \mathbbm{1} \otimes \sigma^z) + h_6 \sigma^z \otimes \sigma^z \nonumber \\
&+ h_7 \sigma^-\otimes \sigma^- + h_8 \sigma^+ \otimes \sigma^+.
\end{align}
can all be mapped to three-site operators.
Note that this is the general form of the ansatz used in Ref.~\cite{de2020classifying}. After a BST, the corresponding three-site operator $\check{G}$ satisfies
\begin{equation}
[\check{G}_{123}, \check{G}_{345}]=0,
\label{comm}
\ee
since the operator acting on the common site (3) is either $\mathbbm{1}$ or $\Sigma^z$.

Since we want to consider open quantum systems, we need to generalize the transformation to mixed density matrices.
The vectorized density matrix, $|\rho\rangle$, is generated by $|s'_1,...\rangle |\tilde{s}'_1,...\rangle\rangle $. In particular, the basis states upon which the two-site gate $\check{R}$ discussed in the previous section acts are generated by $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$.  Since we constructed $\check{R}$ using the Lindbladian of the form in Eq.(\ref{lin}), if both $h$ and $l$ commute with $\sigma^z\otimes \sigma^z$, the action of $\check{R}$ on $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ may preserve the product of $s'_1 s'_2$ and of $\tilde{s}'_1 \tilde{s}'_2$, separately. If this is the case, $|s'_1, s'_2\rangle$ and $ |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ can separately be mapped to $|s_1, s_2, s_3\rangle$ and $|\tilde{s}_1, \tilde{s}_2, \tilde{s}_3\rangle\rangle$. There are four states corresponding to each $|s'_1, s'_2\rangle |\tilde{s}'_1, \tilde{s}'_2\rangle\rangle$ and the (64-dim) Hilbert space separates into four sectors. Moreover, $\check{R}$ can be regarded as a linear operator within each of the sectors. The corresponding three-site gate $\check{G}$ then also preserves each sector. Thus, $\check{G}$ is equal to the direct sum of four copies of $\check{R}$ when acting on the (64-dim) vectorized density operator, which, by linearity, implies that
\begin{align}
\check{R}_{23}(u_1, u_2)&\check{R}_{12}(u_1, u_3) \check{R}_{23}(u_2, u_3) \nonumber \\
&=\check{R}_{12}(u_2, u_3)\check{R}_{23}(u_1, u_3)\check{R}_{12}(u_2, u_3).
\label{RYBE}
\end{align}
leads to
\begin{align}
\check{G}_{234}(u_1, u_2)& \check{G}_{123}(u_1, u_3) \check{G}_{234} (u_2, u_3) \nonumber \\
& = \check{G}_{123}(u_2, u_3) \check{G}_{234}(u_1, u_3) \check{G}_{123} (u_1, u_2).
\label{GYBE}
\end{align}
Also, if two-site (Kraus) operators $K_i$ satisfy
$
K_0^{\dagger}K_0 +K_1^{\dagger}K_1 = \mathbbm{1}
$, the corresponding three-site operators $\tilde{K}_i = K_i\otimes K_i$ satisfy
$
\tilde{K}_0^{\dagger}\tilde{K}_0 +\tilde{K}_1^{\dagger}\tilde{K}_1 = \mathbbm{1}.
$

\begin{figure}[b]
\centering
\includegraphics[width=0.4\textwidth]{circuit3}
\caption{Floquet quantum circuit with staggered special three-site gates. A blue disk represents a general operator while a black dot can be $\mathbbm{1}$ or $\Sigma^z$.}
\label{circuit3}
\end{figure}

\subsection{Integrability structure}
Since the BST is nonlocal, it is  not obvious whether all  local charges of a circuit with two-site gates are mapped to local charges of the corresponding circuit with three-site gates. To show that this is indeed the case, and thus the latter is integrable, we follow Ref. \cite{gombor2021integrable}.

We have already seen that if the two-site gate $\check{R}$ satisfies the YBE (Eq.(\ref{RYBE})) the corresponding three-site gate $\check{G}$ satisfies Eq.(\ref{comm}), Eq.(\ref{GYBE}), $ \check{G}_{123}(v, v)= \mathbbm{1}$, and
\begin{equation} \check{G}_{123}(u, v) \check{G}_{123}(v, u) = \mathbbm{1}.\ee
Define $\check{L}_{123} (u) =\check{G}_{123}(u, 0)$ and
\begin{equation}
\check{R}_{12, 34}(u,v) = \check{L}_{234}(u) \check{G}_{123}(u,v) \check{L}_{234}(v)^{-1}.
\label{rlgl}
\ee
Thus defined four-site $\check{R}$ then satisfies
\begin{align}
&\check{R}_{34, 56}(u_1, u_2) \check{R}_{12, 34}(u_1, u_3) \check{R}_{34, 56} (u_2, u_3) \nonumber \\
& = \check{R}_{12, 34}(u_2, u_3) \check{R}_{34, 56}(u_1, u_3) \check{R}_{12, 34} (u_1, u_2),
\end{align}
and we can construct a transfer matrix
\begin{align}
T(u) &= \text{Tr}_{a,b} [R_{(1,2), (a,b)}(v, u)  R_{(3,4), (a,b)}(v, u)... \nonumber \\
& R_{(L-3,L-2), (a,b)}(v, u) R_{(L-1,L), (a,b)}(v, u)],
\end{align}
where $(a, b)$ denote the ``auxiliary space" and $R_{(a,b), (c,d)} \equiv P_{(a,b), (c,d)}\check{R}_{(a,b), (c,d)}$ ($P$ is the permutation operator).
Regarding each pair of indices as a single index, we can follow the derivation Section \ref{int_str} to show that $ T(v)  = U^{-2}$, where $U$ is the one-site translation operator, and
\begin{equation} [T(u), T(u')] =0.\ee
Moreover, by Eq.(\ref{rlgl}),  $\check{R}_{12, 34}(u,0) = \check{L}_{234}(u) \check{L}_{123}(u)$. Then
\begin{align}
 T(0) &=  \nonumber \\
&\text{Tr}_{a,b} [L_{1,2, b}L_{1,2, a}  L_{3,4, b} L_{3,4, a}\cdots   L_{L-1,L, b}L_{L-1,L, a}],
\end{align}
Define
\begin{equation} {\cal{V}}(v, u) =  T(v)^{-1} T(u),\ee
then the Floquet operator (Fig.\ref{circuit3}) corresponds to
\begin{equation} {\cal{V}}(v, 0) = \prod_{j=1}^{L/2}\check{L}_{2j, 2j+1, 2j+2}(v)\prod_{j=1}^{L/2}\check{L}_{2j-1, 2j, 2j+1}(v).\ee
The operator $ {\cal{V}}(v, u) $  leads to local conserved charges, the first being
\begin{equation}
Q_4  =\frac{d }{d u }  {\cal{V}}(v,u) |_{u =v}.
\ee
It is a sum of four-site operators. Higher derivative yield higher conserved charges. Therefore, we conclude that, if the BST is applicable to an integrable Floquet quantum circuit with two-site gates as discussed above, the corresponding circuit with three-site gates is also integrable. In particular, all five models discussed in the last section can be mapped to integrable circuits with nonunitary three-site gates.

\subsection{Example}
\label{example}
Let us discuss the simplest example:  Model B1 with $\tau \kappa =1$. Both Kraus operators of Model B1 are special cases of a general $\check{R}$-matrix of the zero-field eight-vertex (ZFEV) model:
\begin{equation}
\check{R} =  \begin{pmatrix} a & 0 & 0& d\\ 0 & c & b& 0 \\  0 & b & c& 0 \\ d & 0 & 0& a\end{pmatrix},
\label{R}
\ee
where $a, b, c, d$ are usually functions of spectral parameters $u$ and $v$.  The $\check{R}$-matrix can also be written as
\begin{align}
\check{R} & = \frac{1}{2} (a+c) \mathbbm{1} \otimes \mathbbm{1} +\frac{1}{2} (b+d)  \sigma^x  \otimes  \sigma^x \nonumber \\
&+\frac{1}{2} (b-d)  \sigma^y  \otimes  \sigma^y+ \frac{1}{2}( a- c) \sigma^z \otimes \sigma^z,
\label{R1}
\end{align}
which after a BST yields
\begin{align}
\check{G} & = \frac{1}{2} (a+ d \Sigma_2^x)( 1 + \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} (c+ b \Sigma_2^x)( 1- \Sigma_1^z \Sigma_3^z).
\label{G}
\end{align}
Similarly, the Kraus operators of Model B1 ($\tau \kappa =1$) are mapped to $K_0=  1/\sqrt{u+1} $ and
\begin{equation}
K_1=  \sqrt{\frac{u}{u+1}}[ \frac{1}{2} ( 1 + \Sigma_1^z \Sigma_3^z)+   \frac{1}{2}  \Sigma_2^x( 1- \Sigma_1^z \Sigma_3^z)].
\ee
Their action on the transformed $\rho$ is similar to Eq.(\ref{Kraus}). In the limit $u \to \infty$, the action is unitary with $K_1$ the same as RCA150 (see Appendix \ref{app1}). Other cases in Ref. \cite{de2021constructing} can be described similarly.

\section{Classical SCA}
\label{SCA}
In a standard cellular automaton\cite{Wolfram}, each site is one classical bit and it is updated according to a rule that takes the values of itself and its neighbors as the input \cite{farrelly2020review}. If the update rule is reversible, i.e., each output state has only one predecessor, it is also called a reversible cellular automaton; otherwise, it is irreversible.   Classical SCA can be obtained from QCA by constraining the density matrix $\rho$ to be diagonal, $\rho = \text{diag}(\rho_{\uparrow\uparrow}, \rho_{\uparrow\downarrow}, \rho_{\downarrow\uparrow}, \rho_{\downarrow\downarrow})$. That the restricted $\check{R}$ is stochastic (i.e., it maps a probability distribution to another one) is guaranteed by the CPTP property of the original $\check{R}$. (A unital $\check{R}$ reduces to a bi-stochastic matrix \cite{pavel2021correlations}.) Preservation of the diagonal space leads to classical (discrete) flows of probabilities $P(s_1,s_2,..., s_L) = \langle s_1, s_2,...,s_L| \rho |s_1, s_2,...,s_L\rangle$ and the continuum limit of the five circuits with two-site classical gates is described Ref.\cite{de2021constructing}. As an example, Model B1 with $\tau \kappa =1$, for finite $u$ yields
\begin{align}
\text{diag}&(\rho_{\uparrow\uparrow}, \rho_{\uparrow\downarrow}, \rho_{\downarrow\uparrow}, \rho_{\downarrow\downarrow}) \nonumber \\
&\to \text{diag}(\rho_{\uparrow\uparrow}, \frac{\rho_{\uparrow\downarrow} +u \rho_{\downarrow\uparrow}}{u+1}, \frac{\rho_{\downarrow\uparrow}+u \rho_{\uparrow\downarrow}}{u+1}, \rho_{\downarrow\downarrow}).
\end{align}
The corresponding transition matrix after a BST can be expressed as
\begin{align}
\check{G}(u) & = \frac{1}{2}( 1 + \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} \frac{1+ u \Sigma_2^x}{u+1}( 1- \Sigma_1^z \Sigma_3^z).
\end{align}

These integrable models are a subset of larger integrable family: ZFEV model (cf. Appendix \ref{app2}). Indeed, the application of a BST to the integrable chains Eq.(\ref{R1}) yields Eq.(\ref{G}).
The coefficients are given in Eq.(\ref{abcd}).
For $u \ge 0$ and $v \ge 0$, they are all nonnegative and $a + d = b + c =1$. They correspond to the ZFEV model with
\begin{align}
\Delta  &=  \frac{a^2 + b^2 - c^2 - d^2}{2(a b + c d)} = \text{sech} (2v), \nonumber \\
\Gamma & = \frac{a b - c d}{a b + c d} = \text{sech} (2v).
\end{align}
Since $\text{sech}(2v) \le 1$, the ground state of the ZFEV model is always in the disordered phase or on its critical boundaries \cite{baxter2016exactly}. In particular, RCA150 corresponds to $v=0$ after a BST, i.e., it always sits on a critical line. For a recent analysis, see Ref.\cite{wilkinson2021exact}. For the ZFEV model, there are many symmetries which preserve the integrability \cite{baxter2016exactly}. In particular, $a\leftrightarrow d$ and $b\leftrightarrow c$ maps RCA150 to RCA105. Note that in Ref.\cite{gombor2021integrable}, RCA150 and RCA105 were embedded into two different families of integrable quantum models. This is because the regularity condition was imposed. We show in Appendix \ref{app3} that the regularity condition can be relaxed without breaking the integrability and it corresponds exactly to the $a\leftrightarrow d$ and $b\leftrightarrow c$ symmetry. The transition matrix is mapped to
\begin{align}
\check{G}(u) & = \frac{1}{2}( 1 - \Sigma_1^z \Sigma_3^z)+  \frac{1}{2} \frac{1+ u \Sigma_2^x}{u+1}( 1+ \Sigma_1^z \Sigma_3^z).
\end{align}

Thus, we can see that some integrable SCA can be embedded into a family of integrable nonunitary QCA. In particular, RCA150 (and RCA105) can be embedded into families of integrable SCA, unitary QCA \cite{gombor2021integrable}, and nonunitary QCA.

\section{Conclusion and outlook}
We have constructed several two-site translationally invariant integrable nonunitary quantum circuits using the $\check{R}$-matrices that appeared in constructing integrable Lindblad superoperators in the continuous case  in Ref. \cite{de2021constructing}. The trotterized Hubbard model with imaginary interaction strength \cite{sa2021integrable} is one special case of our construction. We showed that these $\check{R}$-matrices generate integrable open systems by demonstrating that they are CPTP, and thus represent a valid quantum channel. By applying the BST, integrable quantum circuits with three-site gates can also be constructed.  Among them, the classical RCA150 (and RCA105) can be embedded into an integrable family describing quantum open systems. The classical discrete flows on the diagonal of density matrices of the entire family of Model B1 can be embedded into the well-known integrable ZFEV model with parameter restricted such that it ground state is always in the disordered phase (and its boundary critical lines). Since the BST is a special case of the more general Clifford transformations \cite{jones2021integrable}, other descencant integrable nonunitary quantum circuits can be similarly derived.

In general, if we can find CPTP solutions to the YBE, we can use them directly to build integrable nonunitary quantum circuits. No Markovian approximation would be needed. Since we can regard the $\check{R}$-matrices of Model B1 and of Model B2 as those obtained from their unitary analogs via some analytical continuation, a closer look at the analytical continuation of known unitary solutions of YBE may may help find other CPTP solutions.

Our work points to a possible systematic way to obtain CPTP solutions to the YBE. Since all $\check{R}$-matrices of the five models in Ref.\cite{de2021constructing} are CPTP, a natural direction is to verify whether all $\check{R}$-matrices that can be constructed using the Lindbladian approach of Ref. \cite{de2021constructing} satisfy this property. In particular, if we use an ansatz of the form $\check{R}(u) =\mathbbm{1} + \sum_{n =1}^{\infty} r^{(n)}({\cal{L}}) u^n$ where $r^{(n)}({\cal{L}})$ is a polynomial of degree $n$ in ${\cal{L}}$ \cite{de2019classifying}, the trace-preserving property of $\check{R}$ is automatically guaranteed since $\text{Tr}[{\cal{L}}(\cdot)] =0$. Then it is sufficient to show that $\check{R}$ is completely positive.

A natural extension of the result in Ref.\cite{de2021constructing} to multiple families of jump operators may also lead to extension of the results in this paper to more general integrable open quantum circuits. In particular, it would be interesting to know if there is any connection to operator-space fragmentation discussed in Ref. \cite{Essler2020integrability} with two families of jump operators $l_{j, j+1}^{(1)} = \sqrt{J_1} \sigma_j^+\sigma_{j+1}^-$ and $l_{j, j+1}^{(2)} = \sqrt{J_2} \sigma_j^-\sigma_{j+1}^+$.

We have noticed that  both Model B1 and Model B2 have the Kraus representation $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$, a convex combination of unitary $U_i$. It superficially resembles the action of local gates in a random unitary quantum circuit where $U_i$  may be taken to be Haar-random unitary \cite{nahum2017quantum, nahum2018operator, vonkeyserlingk2018operator}  or Clifford-random unitary \cite{li2018quantum, Ippoliti2021post}. Even though random systems are generically chaotic, it is tempting to see if techniques invented in the context of random circuits can be adopted to integrable nonunitary circuits. On the other hand, the action $\phi[\rho] =\sum_i \alpha_i U_i\rho U_i^{\dagger}$ can be interpreted as an ensemble of quantum trajectories \cite{wiseman1996quantum, jian2020criticality}. If we focus on a single quantum trajectory and interpret nonunitary quantum gates as weak measurements following unitary evolution \cite{li2018quantum, Ippoliti2021post}, then we do not have to use vectorized density matrices. In particular, we can construct more integrable non-unitary quantum circuits by applying the spacetime duality to integrable unitary quantum circuits \cite{Ippoliti2021post}. All we need to is relax the inversion property (Eq.(\ref{inv})) and regularity condition (Eq.(\ref{regcon})), which can be easily done. We will leave it for future work.

We have only proved the integrability of the open quantum circuits constructed. How the extensive number of conserved charges dictate the thermalization in integrable nonunitary quantum circuits is an important question. Explicit computation of time correlation functions can be challenging since the gates are no longer unitary. It would be interesting to see if the Bethe ansatz approach can yield useful information about the complex spectrum of the circuits \cite{ziolkowska2020yang} and, moreover, how dominant complex eigenvalues controls the long-time dynamics \cite{claeys2022correlations}. We may draw some inspirations from solutions to integrable systems with boundary drives \cite{vanicat2018integrable}. Numerically study of spectral statistics as in Ref.\cite{sa2021integrable} may provide another perspective. It is also helpful to compute dynamical entanglement-related qualities, such as entanglement spread, entanglement negativity, of the models. It will certainly give us a better understanding of irreversible/nonunitary QCA \cite{farrelly2020review, brennen2003entanglement, piroli2020quantum}.

\section*{Acknowledgement}
We thank Aashish Clerk for discussions. This material is based upon work supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.

redpajama_set_name: RedPajamaArXiv