Qutrit Circuits and Algebraic Relations:
A Pathway to Efficient Spin-1 Hamiltonian Simulation

Oluwadara Ogunkoya
Superconducting Quantum Materials and Systems Center (SQMS), Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
Fermi National Accelerator Laboratory, Batavia, IL, 60510, USA

Joonho Kim
Superconducting Quantum Materials and Systems Center (SQMS), Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
Rigetti Computing, Berkeley, CA 94710, USA

Bo Peng
Physical and Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, WA 99352, USA

A. Barış Özgüler
Superconducting Quantum Materials and Systems Center (SQMS), Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
Fermi National Accelerator Laboratory, Batavia, IL, 60510, USA

Yuri Alexeev
Computational Science Division, Argonne National Laboratory, Lemont, IL 60439, USA

Quantum information processing has witnessed significant advancements through the application of qubit-based techniques within universal gate sets. Recently, exploration beyond the qubit paradigm to $d$-dimensional quantum units or qudits has opened new avenues for improving computational efficiency. This paper delves into the qudit-based approach, particularly addressing the challenges presented in the high-fidelity implementation of qudit-based circuits due to increased complexity. As an innovative approach towards enhancing qudit circuit fidelity, we explore algebraic relations, such as the Yang-Baxter-like turnover equation, that may enable circuit compression and optimization. The paper introduces the turnover relation for the three-qutrit time propagator and its potential use in reducing circuit depth. We further investigate whether this relation can be generalized for higher-dimensional quantum circuits, including a focused study on the one-dimensional spin-1 Heisenberg model. Our work outlines both rigorous and numerically efficient approaches to potentially achieve this generalization, providing a foundation for further explorations in the field of qudit-based quantum computing.

§ INTRODUCTION

Quantum information processing through a gate-based quantum computing approach with qubits involves a universal gate set consisting of single-qubit gates in the $SU(2)$ group and entangling two-qubit gates [21]. This approach has been intensively studied in recent years for applications in quantum information science. For example, in quantum error correction (QEC), although the qubit-based surface code [9, 20] has been thus far the primary route for error detection and correction [17, 40, 37, 54], more favorable QEC schemes are under active exploration.

One alternative to qubit-based quantum information processing is to use $d$-dimensional ($d>2$) quantum units, or qudits. This approach provides an immediate opportunity to explore possibilities beyond two levels as basic units and can therefore boost performance through access to a larger computational space and fewer entangling gates for certain algorithms. The universal gate set for qudits consists of single-qudit gates in the $SU(d)$ group and entangling two-qudit gates [42, 55, 10]. Recently, the universal qudit gate set and their coherent control have been implemented in superconducting transmons [7, 26, 22], photonic circuit [18], and trapped ions [47, 29]. For QEC, qutrits, the simplest form of qudits ($d=3$), have theoretically exhibited advantages through compact logical encoding (to overcome erasure and ternary errors [33, 41, 39]), fault tolerance enhancement [13], and magic state distillation [14]. With regard to other applications of qutrits, there are theoretical reports of utilizing qutrits to improve quantum algorithms (e.g., Shor's factoring [8], Grover's search [25, 12], and quantum Fourier transformation [45]), quantum simulations [28], quantum cryptography [5, 11], and quantum communication [50]. On the other hand, implementing a high-fidelity qudit-based circuits also presents challenges that are essentially associated with increased complexity in the design, fabrication, and control of quantum systems with higher dimensions.

Towards improving the fidelity of qudit circuit, one direction could be to explore some algebraic relations between qudit circuits for the purpose of circuit compression and optimization, which would result in more fault-resilient performance. Similar exploration for the qubit-based circuits have recently become quite active [44, 35, 15, 38, 19, 1, 6, 2]. For example, the Yang-Baxter equation (YBE), which was originally introduced in theoretical physics [52] and statistical mechanics [4], has recently been shown to have connections to topological entanglement, quantum entanglement, and quantum computational universality [24, 43, 34, 53, 51, 3, 46, 27]. In our previous work [46], we proved that for some model systems, the two-qubit time propagator $\mathcal{R}_{\theta,\delta}$, parametrized by a rotation angle $\theta$ and a phase $\delta$, bears a similar algebraic form to the $SU(2)$ solution of the YBE. Therefore, the turn-over relationship (<ref>) can hold as long as certain algebraic relations between the parameters on both sides are satisfied. Remarkably, this turn-over relation can be utilized to compress the corresponding time evolution circuit to a depth that scales linearly with respect to the number of qubits.

\begin{equation}
\resizebox{0.85\hsize}{!}{%
$\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{R}_1(\boldsymbol\alpha)} & \qw & \multigate{1}{\mathcal{R}_3(\boldsymbol\gamma)} & \qw &\\
& \ghost{\mathcal{R}_1(\boldsymbol\alpha)} & \multigate{1}{\mathcal{R}_2(\boldsymbol\beta)} & \ghost{\mathcal{R}_3(\boldsymbol\gamma)} & \qw & ~~~ = \\
& \qw & \ghost{\mathcal{R}_2(\boldsymbol\beta)} & \qw & \qw &}
\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{R}_5(\boldsymbol\epsilon)} & \qw & \qw \\
& \multigate{1}{\mathcal{R}_4(\boldsymbol\delta)} & \ghost{\mathcal{R}_5(\boldsymbol\epsilon)} & \multigate{1}{\mathcal{R}_6(\boldsymbol\zeta)} & \qw \\
& \ghost{\mathcal{R}_4(\boldsymbol\delta)} & \qw & \ghost{\mathcal{R}_6(\boldsymbol\zeta)} & \qw }$
} \label{YBE_general}
\end{equation}

This observation then opens the question of whether these turn-over relations can be generalized for quantum circuits with higher dimensions. Mathematically, there have been discussions on finding the high-dimension solutions to the generalized YBE [49, 23, 48, 16, 30]; however, numerical searching of these high-dimension solutions can be challenging. So far, in addition to the $SU(2)$ solutions, only an $8\times 8$ solution to a generalized YBE has been reported and used to generate the Greenberger-Horne-Zeillinger states [49]. In this paper, as an exploratory effort in this direction, we primarily focus on establishing a similar turn-over relation that can be utilized for the efficient quantum simulation of the quantum time dynamics of the one-dimensional spin-1 Heisenberg model. In particular, we search for (1) the existence of rigorous algebraic conditions for the similar turn-over relations to hold, and (2) a numerically efficient approach that can provide  imprecise but sufficiently accurate turn-over circuits in the absence of rigorous algebraic relations.

In the following sections, we  first define some notations that will be used in this paper. Then, we show for some simple models, rigorous turn-over relations do exist. Finally, for models without rigorous turn-over relations, we propose a numerical recipe to achieve inexact but accurate enough qutrit circuit substitutes. The numerical recipe and the corresponding error analysis are provided for the three-qutrit circuit simulating the time dynamics of a three-site spin-1 Heisenberg model. We conclude this work by offering some remarks on our future effort.

§ NOTATIONS AND SPIN ALGEBRA

The closed-system dynamics of a one-dimensional array of level-$d$ variables is realized by $U(d^N)$ unitary matrices, where $N$ denotes the size of the system. Throughout this paper, however, we specifically treat qudits as spin $s=(d-1)/2$ quantum states and consider their time evolution with certain bilinear spin Hamiltonians.

We recall that for $d=2$, $s=1/2$ the spin operators, satisfying the $SU(2)$ commutation algebra

\begin{align}
[S^x, S^y] = iS^z,~~ [S^y, S^z] = iS^x,~~ [S^z,S^x] = iS^y,
\end{align}

are halves of the Pauli matrices:

\begin{align}
X = \left( \begin{array}{cc} 0 &  1 \\ 1 &  0 \end{array} \right),~~
Y = \left( \begin{array}{cc} 0 & -i \\ i &  0 \end{array} \right),~~
Z = \left( \begin{array}{cc} 1 &  0 \\ 0 & -1 \end{array} \right).
\end{align}

The $SU(2)$ algebra allows a  quadratic Casimir invariant,

\begin{align}
(S^x)^2 + (S^y)^2 + (S^z)^2 =  s(s+1) \mathbf{1}. %= \tfrac{d^2-1}{4} \mathbf{1}.
\end{align}

Under the $z$-basis representation of qudit states, i.e.,

\begin{align}
\label{eq:spin-z}
S^z|l\rangle = (l-(d-1)/2)|l\rangle \quad \text{with}\quad 0 \leq l < d,
\end{align}

the $d\times d$ matrix representation of $S^x$ and $S^y$ are

\begin{align}
\label{eq:spin-x}
(S^x)_{pq} &=
\begin{dcases}
\tfrac{1}{2} \lambda_{\text{min(p, q)}} & \text{if $|p-q| = 1$}\\
0 & \text{otherwise}
\end{dcases}\\
\label{eq:spin-y}
(S^y)_{pq} &=
\begin{dcases}
+\tfrac{1}{2i} \lambda_{\text{min(p, q)}} & \text{if $p-q = +1$}\\
-\tfrac{1}{2i} \lambda_{\text{min(p, q)}} & \text{if $p-q = -1$}\\
0 & \text{otherwise}
\end{dcases}
\end{align}

with $\lambda_l \equiv \sqrt{(d-l)l}$. When restricted to binary systems, they are the Pauli matrices satisfying the Clifford algebra $\{S^a, S^b\} = 2 \delta^{ab}$. Note that the same anti-commutation relation no longer holds for $d>2$, thus, e.g., $(S^a)^2 \neq \mathbf{1}.$

For three-level systems ($d=3$, $s=1$) the $z$-basis representation of the spin-$1$ operators becomes

\begin{align}
S^x &= \frac{1}{\sqrt{2}}\begin{pmatrix}0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0\end{pmatrix}
= \frac{1}{\sqrt{2}}\left( X \oplus 0 + 0 \oplus X\right), \label{eq:qutrit_x} \\
S^y &= \frac{1}{\sqrt{2}}\begin{pmatrix}0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0\end{pmatrix}
= \frac{1}{\sqrt{2}}\left( Y \oplus 0 + 0 \oplus Y\right), \label{eq:qutrit_y} \\
S^z &= \begin{pmatrix}1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1\end{pmatrix}
= \left( Z \oplus 0 + 0 \oplus Z\right). \label{eq:qutrit_z}
\end{align}

As an alternative to (<ref>), (<ref>), and (<ref>), it is sometimes more convenient to use the adjoint representation of the spin operators,

\begin{align}
\tilde{S}^x &= \begin{pmatrix}0 & 0 & 0 \\  0 & 0 & i \\  0 & -i & 0\end{pmatrix} = 0 \oplus (-Y), \label{sx_adj} \\
\tilde{S}^y &= \begin{pmatrix}0 & 0 & i \\  0 & 0 & 0 \\ -i &  0 & 0\end{pmatrix}
= P_y \tilde{S}^x P_y^\dagger, \label{sy_adj} \\
\tilde{S}^z &= \begin{pmatrix}0 & i & 0 \\ -i & 0 & 0 \\  0 &  0 & 0\end{pmatrix}
= P_z \tilde{S}^x P_z^\dagger \label{sz_adj}
\end{align}

with permutation matrices $P_y$ and $P_z$ given by

\begin{align}
P_y = \begin{pmatrix}0&1&0\\ 1&0&0\\ 0&0&1 \end{pmatrix},~~
P_z = \begin{pmatrix}0&1&0\\ 0&0&1\\ 1&0&0 \end{pmatrix},
\end{align}

\begin{align}
(\tilde{S}^x)^2 &= \text{diag}(0,1,1),\\
% \begin{pmatrix}0 & 0 & 0 \\  0 & 1 & 0 \\  0 & 0 & 1\end{pmatrix},\\
(\tilde{S}^y)^2 &=
\text{diag}(1,0,1),\\
% \begin{pmatrix}1 & 0 & 0 \\  0 & 0 & 0 \\  0 & 0 & 1\end{pmatrix},\\
(\tilde{S}^z)^2 &=
\text{diag}(1,1,0).
% \begin{pmatrix}1 & 0 & 0 \\  0 & 1 & 0 \\  0 & 0 & 0 \end{pmatrix},
\end{align}

Similar to the $SU(2)$ cases, $\{\tilde{S}^x,\tilde{S}^y,\tilde{S}^z\}$ follows the commutation algebra

\begin{align}
[\tilde{S}^x, \tilde{S}^y] = i\tilde{S}^z,~~ [\tilde{S}^y, \tilde{S}^z] = i\tilde{S}^x,~~ [\tilde{S}^z, \tilde{S}^x] = i\tilde{S}^y
\end{align}

The two representations $\{S^x,S^y,S^z\}$ and $\{\tilde{S}^x,\tilde{S}^y,\tilde{S}^z\}$ are connected through the basis change between spherical and Cartesian coordinates,

\begin{align}
\label{eq:qutrit-basis-change}
\tilde{S}^a = U_{+} S^a U^\dagger_{+},~~ a = x, y, z
\end{align}

\begin{align}
U_\pm = \frac{1}{\sqrt{2}}\begin{pmatrix}
-1 & 0 & \pm 1 \\
+i & 0 & \pm i \\
0 & \sqrt{2} & 0
\end{pmatrix}, \label{mat_U}
\end{align}

The basis change (<ref>) does not affect the algebraic relations and the circuit substitutes established in the following sections.

Some algebraic relations of the spin-1 operators are worth mentioning. For example, for $n\ge 1$, we have

\begin{align}
\begin{cases}
(\tilde{S}^a)^{2n} = (\tilde{S}^a)^2\\ (\tilde{S}^a)^{2n+1} = \tilde{S}^a \\ \tilde{S}^a\tilde{S}^b\tilde{S}^a = \mathbf{0}_3
\end{cases} \quad \text{for } a\neq b\in \{x,y,z\}.
\end{align}

which implies

\begin{align}
\mathcal{U}_x(\alpha) &= \exp{(-i \alpha \tilde{S}^x \otimes \tilde{S}^x)} \label{eq:Ux} \\
&= \mathbf{I}_9 - i\sin(\alpha)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\alpha}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2  \notag \\
&= \begin{pmatrix}
\mathbf{I}_3 & \mathbf{0}_3 & \mathbf{0}_3 \\
\mathbf{0}_3 & \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^x)^2 & \sin(\alpha)\tilde{S}^x \\
\mathbf{0}_3 & -\sin(\alpha)\tilde{S}^x & \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^x)^2
\end{pmatrix}\notag,
\end{align}

\begin{align}
\mathcal{U}_y(\alpha) &= \exp{(-i \alpha \tilde{S}^y \otimes \tilde{S}^y)} \label{eq:Uy} \\
&= \mathbf{I}_9 - i\sin(\alpha)(\tilde{S}^y \otimes \tilde{S}^y) - 2\sin^2(\frac{\alpha}{2})(\tilde{S}^y \otimes \tilde{S}^y)^2  \notag \\
&= \begin{pmatrix}
\mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^y)^2 & \mathbf{0}_3 & \sin(\alpha)\tilde{S}^y \\
\mathbf{0}_3 & \mathbf{I}_3 & \mathbf{0}_3 \\
-\sin(\alpha)\tilde{S}^y & \mathbf{0}_3 & \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^y)^2
\end{pmatrix}\notag,
\end{align}

\begin{align}
\mathcal{U}_z(\alpha) &= \exp{(-i \alpha \tilde{S}^z \otimes \tilde{S}^z)} \label{eq:Uz} \\
&= \mathbf{I}_9 - i\sin(\alpha)(\tilde{S}^z \otimes \tilde{S}^z) - 2\sin^2(\frac{\alpha}{2})(\tilde{S}^z \otimes \tilde{S}^z)^2  \notag \\
&= \begin{pmatrix}
\mathbf{I}_3 - 2(\sin(\frac{\alpha}{2})\tilde{S}^z)^2 & \sin(\alpha)\tilde{S}^z & \mathbf{0}_3 \\
-\sin(\alpha)\tilde{S}^z & \mathbf{I}_3 - 2(\sin(\frac{\alpha}{2})\tilde{S}^z)^2 & \mathbf{0}_3 \\
\mathbf{0}_3 & \mathbf{0}_3 & \mathbf{I}_3
\end{pmatrix}, \notag
\end{align}

where $\mathbf{I}_m$ denotes an $m\times m$ identity matrix and $\mathbf{0}_m$ denotes an $m\times m$ zero matrix. Notably, from (<ref>) and (<ref>)

\begin{align}
&\Qcircuit @C=.5em @R=.5em {
& \multigate{2}{\mathcal{U}_y(\alpha)} & \qw &\\
& & & ~~ = \\
& \ghost{\mathcal{U}_y(\alpha)} & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \gate{P_y} & \multigate{1}{\mathcal{U}_x(\alpha)} & \gate{P_y{{{}^\dagger}}} & \qw \\
& \gate{P_y} & \ghost{\mathcal{U}_x(\alpha)} & \gate{P_y{{{}^\dagger}}} & \qw } \label{Uy} \\
&\Qcircuit @C=.5em @R=.5em {
& \multigate{2}{\mathcal{U}_z(\alpha)} & \qw & \\
& & & ~~ = \\
& \ghost{\mathcal{U}_z(\alpha)} & \qw & }
~~~~\Qcircuit @C=.5em @R=.5em {
& \gate{P_z} & \multigate{1}{\mathcal{U}_x(\alpha)} & \gate{P_z{{{}^\dagger}}} & \qw \\
& \gate{P_z} & \ghost{\mathcal{U}_x(\alpha)} & \gate{P_z{{{}^\dagger}}} & \qw } \label{Uz}
\end{align}

A more interesting feature of $\mathcal{U}_{a}(\alpha)$ with $a\in\{x,y,z\}$ is that

\begin{align}
&\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & \qw &\\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & ~~= \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw \\
& \multigate{1}{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw \\
& \ghost{\mathcal{U}_a(\alpha)} & \qw & \qw } \label{turnover1}
\end{align}

§ YANG-BAXTER-LIKE RELATIONS IN QUTRIT CIRCUIT

This section aims to search Yang-Baxter-type identities for qutrit circuits. Specifically, this means establishing a $(3^3 \times 3^3)$ matrix relation of the following type:

\begin{align}
\begin{split}
&\big(\mathcal{R}_1(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{R}_2(\beta)\big)\big(\mathcal{R}_3(\gamma) \otimes \mathbf{I}_3\big) \\ &~~~~~~~~=\big(\mathbf{I}_3 \otimes \mathcal{R}_4(\delta)\big)\big(\mathcal{R}_5(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{R}_6(\zeta)\big)
\end{split}
\label{eq:ybe-qutrit}
\end{align}

where the Greek letters $\alpha, \cdots, \zeta$ collectively denote continuous rotations that parameterize two-qutrit operators $\mathcal{R}_n$, defined as a product of $\mathcal{U}_a$'s ($a\in \{x, y, z\}$).

The rotation angles on the LHS of (<ref>) are unrestricted; we require that for all values of $\alpha$, $\beta$, $\gamma$, there should be a value of $\delta$, $\epsilon$, $\zeta$ that satisfies (<ref>).
The parameters are typically related via unclosed expressions involving trigonometric functions, derived from element-wise equalities of (<ref>).

In the following subsections, we describe analytical and numerical methods to establish (<ref>) with various $\mathcal{R}_n$'s.

§.§ Simple turn-over identities

We start from the simple case where $\mathcal{R}_n = \mathcal{U}_a$ for $a\in \{x, y, z\}$. If all the rotations are same, then by directly applying (<ref>) to the LHS or RHS of (<ref>), we obtained the following YBEs relationships:

$\bullet$ LHS:

\begin{align}
\resizebox{\hsize}{!}{%
$\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw &\\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw & ~~= \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{U}_a(2\alpha)} & \qw & \qw &\\
& \ghost{\mathcal{U}_a(2\alpha)} & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & ~~= \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{U}_a(2\alpha)} & \qw \\
& \multigate{1}{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(2\alpha)} & \qw \\
& \ghost{\mathcal{U}_a(\alpha)} & \qw & \qw }$
} \notag
\end{align}

$\bullet$ RHS:

\begin{align}
\resizebox{\hsize}{!}{%
$\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & \qw &\\
& \multigate{1}{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & ~~= \\
& \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & \qw &\\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_a(2\alpha)} & \qw & ~~= \\
& \qw & \ghost{\mathcal{U}_a(2\alpha)} & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw \\
& \multigate{1}{\mathcal{U}_a(2\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw \\
& \ghost{\mathcal{U}_a(2\alpha)} & \qw & \qw }$
} \notag
\end{align}

$\bullet$ LHS $=$ RHS:

\begin{align}
\resizebox{\hsize}{!}{%
$\Qcircuit @C=.5em @R=.5em {
& \multigate{1}{\mathcal{U}_a(\alpha)} & \qw & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw &\\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_a(2\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw & ~~= \\
& \qw & \ghost{\mathcal{U}_a(2\alpha)} & \qw & \qw &}
~~~~\Qcircuit @C=.5em @R=.5em {
& \qw & \multigate{1}{\mathcal{U}_a(2\alpha)} & \qw & \qw \\
& \multigate{1}{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(2\alpha)} & \multigate{1}{\mathcal{U}_a(\alpha)} & \qw  \\
& \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw }$
} \notag
\end{align}

For more general cases where the rotations are not necessarily same, as explicitly shown in Appendix <ref>, as long as the following conditions are satisfied

\begin{align}
\label{eq:param}
\bigg\{\begin{array}{l}
\alpha + \gamma = \epsilon + 2k\pi  \\
\delta + \zeta = \beta + 2k\pi
\end{array},
~~ k\in \mathbf{Z},
\end{align}

the identity (<ref>) holds for all $\mathcal{U}_a$ $(a = x, y, z)$.
This can also be understood as an immediate consequence of the qubit relations. See Appendix <ref> for the explanation.

\begin{align}
\label{eq:simple-ybe-additional-relation}
\begin{split}
s_{\beta} &= s_{\delta+\zeta}, \quad s_{\epsilon} = s_{\alpha+\gamma}\\
c_{\beta} &= s_{\delta+\zeta}, \quad c_{\epsilon} = c_{\alpha+\gamma}.
\end{split}
\end{align}

Recall that the angle constraints on the qubit circuit have the two solutions <cit.>,

\begin{align}
\begin{dcases}
\beta = \delta+\zeta+\pi\\
\epsilon = \alpha+\gamma+\pi
\end{dcases}, \quad \begin{dcases}
\beta = \delta+\zeta\\
\epsilon = \alpha+\gamma
\end{dcases}
\end{align}

Only the second solution satisfies (<ref>). This shows that

The bonus identities can also be obtained from a simple observation: The spin matrices expressed in the adjoint representation also satisfy, for example,

\begin{align}\label{eq:Sy}
U^\dagger_\pm \tilde{S}^y U_\pm &= \tfrac{1}{\sqrt{2}}(-\tilde{S}^z \mp \tilde{S}^x),
\end{align}

Other cyclic relations can further be obtained by utilizing (<ref>) and (<ref>).
It can be shown that

\begin{align}
V^\dagger_\pm \tilde{S}^x V_\pm &= \tfrac{1}{\sqrt{2}}(-\tilde{S}^y \pm \tilde{S}^z),\label{eq:Sx}\\
W^\dagger_\pm \tilde{S}^z W_\pm &= \tfrac{1}{\sqrt{2}}(+\tilde{S}^x \mp \tilde{S}^y).
\end{align}

\begin{align}
V_\pm &= \frac{1}{\sqrt{2}}\begin{pmatrix}
0 & +i & \pm i \\
0 & \mp 1 & + 1 \\
\sqrt{2} & 0 & 0 \\
\end{pmatrix}, \label{mat_V}\\
W_\pm &= \frac{1}{\sqrt{2}}\begin{pmatrix}
\mp 1 & + 1 & 0 \\
0 & 0 & \sqrt{2} \\
+i & \pm i & 0 \\
\end{pmatrix}. \label{mat_W}
\end{align}

Since the Yang-Baxter identity (<ref>) is independent of the basis, it implies that $\mathcal{R}_n(\theta) \equiv \exp(i\theta (S^a \pm S^b) \otimes (S^a \pm S^b))$ with any $a \neq b \in \{x,y,z\}$ also satisfies (<ref>).

In the same way, we can derive the identity (<ref>) for the following set of extended operators: $\mathcal{R}_n(\theta) = \exp(i\theta (S^x \pm S^y \pm S^z) \otimes (S^x \pm S^y \pm S^z)))$ with any choice of $\pm$ factors. This starts with observing the conjugation relation,

\begin{align}
M^\dagger_{\pm\pm}\, \tilde{S}^z \, M_{\pm\pm} &= \tfrac{1}{\sqrt 3}(\tilde{S}^x \mp \tilde{S}^y \pm \tilde{S}^z),
\end{align}

\begin{align}
\label{eq:conj-matrix-3}
M_{s_1s_2} = \frac{1}{\sqrt{6}}\begin{pmatrix}
-s_1\, \sqrt{3}i & \sqrt{3} i & 0 \\
-s_2\, i  & -s_1 s_2 \,i & 2 i \\
s_2\sqrt{2} & s_1s_2\sqrt{2} & \sqrt{2} \\
\end{pmatrix}
\end{align}

for $s_1, s_2 \in \{\pm 1\}$. Note that the unitarity holds for $M_{\pm\pm}$,
as for the other matrices in (<ref>), (<ref>), and (<ref>).
Thus, any fixed combination of $(S^x \pm S^y \pm S^z)$  can replace all $S^z$ that appear in the circuit identity (<ref>) through $\mathcal{R}_n = \mathcal{U}_z$ if we insert the resolution of identity

\begin{align*}
M_{\pm\pm}M_{\pm\pm}^\dagger = \mathbf{I}_3
\end{align*}

wherever needed.

The discussion so far has established the Yang-Baxter-type circuit relations (<ref>) for the simple cases
where it is assumed that all $\mathcal{R}$ operators have the same form

\begin{align}
\mathcal{R}_n = \mathcal{U}_{s_xs_ys_z}
\end{align}

and depend on a single continuous parameter, i.e.,

\begin{gather}
\mathcal{U}_{s_xs_ys_z} (\alpha) \equiv \exp{\Big(-i\alpha(\textstyle\sum_{a \in \{x,y,z\}} s_a S^a)^{\otimes 2} \Big)}
\end{gather}

where $s_x, s_y, s_z \in \{+1, 0, -1\}$.

§.§ Numerical methods for approximate identities

Unitary Pairs on Test ($W_L \rightarrow W_R$)

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_x(\theta)} & \qw & \multigate{1}{\mathcal{U}_y(\theta)}& \qw & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_x(\theta)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)}& \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \qw & \multigate{1}{\mathcal{U}_y(\zeta)}& \qw & \multigate{1}{\mathcal{U}_x(\mu)} & \qw\\
& \multigate{1}{\mathcal{U}_y(\lambda)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_x(\sigma)} & \ghost{\mathcal{U}_a(\alpha)} & \qw\\
&  \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)}& \qw & \qw\\ } \\
\end{align*}

Trotter scheme $T_1$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_x(\theta)} & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \qw & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_x(\theta)} & \multigate{1}{\mathcal{U}_y(\theta)}  & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \qw & \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)}& \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \qw & \qw& \multigate{1}{\mathcal{U}_y(\zeta)} & \multigate{1}{\mathcal{U}_x(\mu)} & \qw\\
& \multigate{1}{\mathcal{U}_y(\lambda)} & \multigate{1}{\mathcal{U}_x(\sigma)} & \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw\\
&  \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)} & \qw & \qw & \qw\\ } \\
\end{align*}

Trotter scheme $T_2$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_{x+y}(\theta)} & \qw & \qw & \\
& \ghost{\mathcal{U}_{x+y}(\alpha)} & \multigate{1}{\mathcal{U}_{x+y}(\theta)}  & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_{x+y}(\alpha)} & \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \qw& \multigate{1}{\mathcal{U}_{x+y}(\mu)} & \qw\\
& \multigate{1}{\mathcal{U}_{x+y}(\lambda)} & \ghost{\mathcal{U}_{x+y}(\alpha)} & \qw\\
& \ghost{\mathcal{U}_{x+y}(\alpha)} & \qw & \qw\\ } \\
\end{align*}

Trotter scheme $T_3$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_x(\theta)} & \qw & \qw & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_y(\theta)} & \multigate{1}{\mathcal{U}_x(\theta)} & \ghost{\mathcal{U}_a(\alpha)} & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)}& \qw & \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_y(\lambda)} & \qw & \qw & \multigate{1}{\mathcal{U}_x(\mu)} & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_x(\sigma)} & \multigate{1}{\mathcal{U}_y(\zeta)} & \ghost{\mathcal{U}_a(\alpha)} & \qw &\\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)}& \qw & \qw & \\ }  \\
\end{align*}

Trotter scheme $T_4$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_{x+y}(\theta)} & \qw & \qw & \qw & \\
& \ghost{\mathcal{U}_{a+a}(\alpha)} & \multigate{1}{\mathcal{U}_x(\theta)} & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)}&  \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \qw & \qw & \multigate{1}{\mathcal{U}_{x+y}(\mu)} & \qw & \\
& \multigate{1}{\mathcal{U}_{y}(\lambda)} & \multigate{1}{\mathcal{U}_x(\sigma)} & \ghost{\mathcal{U}_{a+a}(\alpha)} & \qw &\\
& \ghost{\mathcal{U}_a(\alpha)} & \ghost{\mathcal{U}_a(\alpha)}& \qw & \qw & \\ }  \\
\end{align*}

Trotter scheme $T_5$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_x(\theta)} & \qw  & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_{x+y}(\theta)}  & \ghost{\mathcal{U}_a(\alpha)} & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_{a+a}(\alpha)} &  \qw & \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_y(\lambda)} & \qw & \multigate{1}{\mathcal{U}_x(\mu)} & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_{x+y}(\mu)}  & \ghost{\mathcal{U}_a(\alpha)} & \qw &\\
& \qw & \ghost{\mathcal{U}_{a+a}(\alpha)} & \qw & \qw & \\ }  \\
\end{align*}

Trotter scheme $T_6$

\begin{align*}
\Qcircuit @C=.5em @R=.5em{
& \multigate{1}{\mathcal{U}_x(\theta)} & \qw & \multigate{1}{\mathcal{U}_y(\theta)}& \qw & \qw & \\
& \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_x(\theta)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_y(\theta)} & \qw & \hspace{0.5cm} \longrightarrow \\
& \qw & \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)}& \qw & \\ }
~~\hspace{0.5cm}\Qcircuit @C=.5em @R=.5em{
& \qw & \multigate{1}{\mathcal{U}_x(\zeta)}& \qw & \multigate{1}{\mathcal{U}_y(\mu)} & \qw\\
& \multigate{1}{\mathcal{U}_x(\lambda)} & \ghost{\mathcal{U}_a(\alpha)} & \multigate{1}{\mathcal{U}_y(\sigma)} & \ghost{\mathcal{U}_a(\alpha)} & \qw\\
&  \ghost{\mathcal{U}_a(\alpha)} & \qw & \ghost{\mathcal{U}_a(\alpha)}& \qw & \qw\\ } \\
\end{align*}

An exact-in-principle identity used for setting up the lower bounds in the numerical tests.

A list of qutrit unitary pairs $W_L(\theta, \cdots, \theta)$ and $W_R(\bm{\theta}_R)$, for which we test the approximate circuit relations through numerical minimization of the infidelity function (<ref>). The top 6 rows correspond to the circuit reflection for different Trotter schemes (<ref>)-(<ref>). While the unitaries in the bottom row are shown to be identical through (<ref>), under the parameter constraint of (<ref>), their equivalence is also tested numerically as a proxy to measure numerical deviations and limitations of optimizers.

We remark that when establishing the circuit identity (<ref>) analytically, numerical validation serves as a useful and efficient strategy to verify its correctness. This step involves evaluating the following expression:

\begin{align}
\begin{split}
\Vert &\big(\mathcal{R}_1(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{R}_2(\beta)\big)\big(\mathcal{R}_3(\gamma) \otimes \mathbf{I}_3\big) \\&-\big(\mathbf{I}_3 \otimes \mathcal{R}_4(\delta)\big)\big(\mathcal{R}_5(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{R}_6(\zeta)\big) \Vert < \varepsilon
\end{split}
\end{align}

with a sufficiently small value of $\varepsilon$. The values for $\alpha$, $\beta$, $\gamma$, and $\delta$ are repeatedly sampled from a uniform distribution while $\epsilon$ and $\zeta$ follow from the constraint (<ref>).

This numerical approach holds a broader range of applications and offers advantages in establishing approximate identities. While exact turnover relations may be specific to certain Hamiltonians, there are chances to develop approximate relations with tolerable levels of infidelity for a wider class of Hamiltonians. Such relations can lead to the compression of circuit depth, thus enhancing the overall fidelity of Trotter circuits running on imperfect hardware. We explore this scenario in the current section, using a concrete example of the spin-1 $XY$ Hamiltonian on three qutrits.

\begin{align}
H_{\rm XY} = -J\sum_{i=0}^2 \left(\tilde{S}_i^x \tilde{S}_{i+1}^x + \tilde{S}_i^y \tilde{S}_{i+1}^y\right).
\label{eq:Ham}
\end{align}

The time evolution unitary of this Hamiltonian system can be written as

\begin{align}
e^{-i t H_{\rm XY}}
&= e^{i Jt \sum_{i=0}^2 (\tilde{S}_i^x \otimes \tilde{S}_{i+1}^x + \tilde{S}_i^y \otimes \tilde{S}_{i+1}^y)}
\label{eq:H_XY}
\end{align}

for which we consider a few available Trotter forms, and numerically explore if the approximate turn-over relation (<ref>) holds. Specifically, for each Trotterization scheme, we evaluate $W_L(\bm{\theta}_L)$ and $W_R(\bm{\theta}_R)$, both representing unitary operators for one Trotter step and its corresponding turnover counterpart, respectively. Afterwards, we minimize the infidelity between $W_L$ and $W_R$,

\begin{align}
\label{eq:CostFunction}
\mathcal{C}(\bm{\theta}_L, \bm{\theta}_R)= 1- \frac{1}{(3^3)^2}\left\Vert\text{tr}\,{(W_L W_R^\dagger)}\right\Vert^2,
\end{align}

by optimizing $\bm{\theta}_R$ for randomly selected values of $\bm{\theta}_L$.

Note that the error analysis in (<ref>) is state-independent, and the fidelity function $\frac{1}{3^6}\Vert\text{tr}\, {(W_L W_R^\dagger)}\Vert^2$ is equivalent to the mean overlap function which is a reasonable figure of merit to quantify coherent errors.

The specific configurations of unitary pairs, $(W_L, W_R)$, that we compute to minimize the infidelity (<ref>) for, are summarized in Fig. <ref>. Each row in the figure is related to the respective Trotter scheme for (<ref>), described below:

\begin{align}
% \begin{split}
T_1 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{y}^{1,2}(\theta) \,
\mathcal{U}_{y}^{0,1} (\theta) \,
\mathcal{U}_{x}^{1,2}(\theta) \,
\mathcal{U}_{x}^{0,1} (\theta)
\Big)^n\label{eq:trot1}\\
T_2 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{y}^{1,2}(\theta) \,
\mathcal{U}_{x}^{1,2}(\theta) \,
\mathcal{U}_{y}^{0,1} (\theta) \,
\mathcal{U}_{x}^{0,1} (\theta)
\Big)^n \label{eq:trot2} \\
T_3 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{x+y}^{1,2}(\theta)\,
\mathcal{U}_{x+y}^{0,1}(\theta)
\Big)^n\label{eq:trot3}\\
T_4 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{y}^{0,1}(\theta)\,
\mathcal{U}_{x}^{1,2}(\theta)\,
\mathcal{U}_{y}^{1,2}(\theta)\,
\mathcal{U}_{x}^{0,1}(\theta)
\Big)^n\label{eq:trot4}\\
T_5 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{y}^{1,2}(\theta)\,
\mathcal{U}_{x}^{1,2}(\theta)\,
\mathcal{U}_{x+y}^{0,1}(\theta)
\Big)^n\label{eq:trot5}\\
T_6 &= \lim_{n\rightarrow \infty} \Big(
\mathcal{U}_{y}^{0,1}(\theta)\,
\mathcal{U}_{x+y}^{1,2}(\theta)\,
\mathcal{U}_{x}^{0,1}(\theta)
\Big)^n \label{eq:trot6}
% \end{split}
\end{align}

where $\theta = t / n$ and

\begin{align}
\begin{split}
\mathcal{U}_{x}^{i,j}(\theta) &= \exp\big({-i\theta \tilde{S}^x_i\otimes\tilde{S}^x_j}\big)\\
\mathcal{U}_{y}^{i,j}(\theta) &= \exp\big({-i\theta \tilde{S}^y_i\otimes\tilde{S}^y_j}\big)\\
\mathcal{U}_{x+y}^{i,j}(\theta) &= \exp\big({-i\theta (\tilde{S}^x_i\otimes\tilde{S}^x_j + \tilde{S}^y_i\otimes\tilde{S}^y_j)}\big).
\end{split}
\end{align}

Although the Trotter unitaries live on a one-dimensional slice of the $\bm{\theta}_L$-parameter space, i.e., $\bm{\theta}_L = (\theta, \theta, \cdots)$, it is necessary to treat all components of $\bm{\theta}_R$ independently in order to achieve a reasonably high fidelity. Therefore, we conduct numerical optimization for $\bm{\theta}_R$ with the following configuration:

\begin{align*}
\min_{\bm{\theta}_R}\mathcal{C}(\theta, \cdots, \theta, \bm{\theta}_R) = \max_{\bm{\theta}_R} \Vert\text{tr}\,{[W_L(\theta, \cdots, \theta)\,W_R^\dagger(\bm{\theta}_R)]}\Vert^2.
\end{align*}

More generally, we consider numerical optimization using respective products of $n_b>1$ instances of $W_L$ and $W_R$.

The minimized infidelity, $\log_{10}\min_{\bm{\theta}_R}\mathcal{C}(\theta, \cdots, \theta, \bm{\theta}_R)$, is obtained through parameter optimization of $\bm{\theta}_R$ across various spin couplings $J \in \{0.1, 0.55, 1.0\}$ and Trotter schemes $T \in \{T_1, T_2, \cdots, T_6\}$. The displayed values are on the logarithmic base 10 scale. Circuit diagrams of parameterized unitaries $W_L(\theta, \cdots, \theta)$ and $W_R(\bm{\theta}_R)$ for each Trotter scheme are shown in Fig. <ref>. The parameter optimization was performed using the BFGS algorithm. We consider the minimized infidelity to be reasonably low if it closely matches the `lower bound', which solely accounts for the numerical inaccuracy of the exact identity (<ref>).

Our optimization results for infidelity are depicted in Fig. <ref> across different values of spin-spin coupling $J \in \{0.10, 0.55, 1.00\}$. Each panel is associated with specific $J$ values and displays the minimized infidelities on a logarithmic scale for different Trotter schemes $\{T_1, \cdots, T_6\}$, covering a range of repetition numbers $1 \leq n_b \leq 5$.
When benchmarking the infidelities between mirror-symmetric pairs of different candidate unitaries, it is necessary to establish a lower bound result that can represent practical expectations for the optimum, taking into account numerical deviations and limitations of optimizers. We achieve this through the following exact-in-principle circuit identity,

\begin{align}
\mathcal{U}_{y}^{0,1} (\theta) \,
\mathcal{U}_{x}^{1,2}(\theta) \,
\mathcal{U}_{x}^{0,1} (\theta) \notag \\
& ~~~~~~~~ =
\mathcal{U}_{y}^{0,1} (\mu) \,
\mathcal{U}_{y}^{1,2}(\sigma) \,
\mathcal{U}_{x}^{0,1} (\zeta)\,
\mathcal{U}_{x}^{1,2}(\lambda), \label{eq:lbound}
\end{align}

where $\mu$, $\sigma$, $\zeta$, and $\lambda$ should in principle be equal to $\theta$,
illustrated in the bottom row of Fig. <ref>. It is derived from the repeated application of the exact turn-over relation (<ref>). Therefore, its minimized infidelities should vanish ideally, i.e., $\min_{\bm{\theta}_R}\mathcal{C}(\bm{\theta}_L, \bm{\theta}_R) = 0$ for any $\bm{\theta}_L$, but are realistically sustained at values ranging from $10^{-7}$ to $10^{-14}$ across different setups of $J$ and $n_b$, as depicted in Fig. <ref>.

In the upper six rows of Fig. <ref>, the reflection pairs of unitaries are strategically arranged so that their repeated application within the Trotter circuit leads to a substantial reduction in the total gate count.  For instance, let us consider the Trotter scheme $T_1$ and its corresponding unitaries. The initial configuration involves 4 $\mathcal{U}_x$ and $\mathcal{U}_y$ operations for every Trotter step, summing up to a total of $4n$ gates. However, by replacing every alternating $W_L(\theta, \cdots, \theta)$ with $W_R(\mu, \sigma, \zeta, \lambda)$ at precomputed values of $\mu$, $\sigma$, $\zeta$, and $\lambda$, it becomes possible to condense $(n-1)$ gates within intermediate Trotter layers due to subsequent applications of the same unitaries, namely,

\begin{align}
\mathcal{U}_a^{i,j}(\theta_1)\mathcal{U}_a^{i,j}(\theta_2) = \mathcal{U}_a^{i,j}(\theta_1 + \theta_2) \text{ with } a \in \{x, y\}.
\end{align}

This results in a reduced count of $(3n+1)$ gates.

Hence, to attain a computational advantage with the Trotter unitary, we select $T \in \{T_1, \cdots, T_6\}$ and $1 \leq n_b \leq 5$ based on initial calculations in Fig. <ref>, then apply the approximate relation $W_L^{n_b} \simeq W_R^{n_b}$ for every alternate set of $n_b$ Trotter steps. While replacing $W_L^{n_b}$ with $W_R^{n_b}$ may increase the overall infidelity of the Trotter circuit, the approximation error has a negligible impact when we use $T$ and $n_b$ whose corresponding infidelity from Fig. <ref> stays at a level similar to the lower bound (<ref>).

The return probability (<ref>) of the spin-1 XY model on a three-site lattice, starting and ending at the state $|202\rangle$, is shown as a function of time within $0 \leq t < 5$. The spin coupling is set at $J=0.55$. The dynamic simulation is represented as the blue line. Additionally, we employed Trotterization of the time-evolution operator (<ref>) over 200 steps, with a corresponding step size of $\theta = 0.025$. Specifically, the Trotter scheme $T_3$ is used, and the resulting data points are indicated by small red dots.
Finally, we apply the circuit compression strategy detailed in Section <ref> to reduce 132 $\mathcal{U}_{x+y}$ gates. More generally,
for a Trotter circuit with $n>3$ steps, the approximate count of reduced gates is $2n/3$. The outcomes from the approximately compressed circuits are denoted by green triangles.

As a pilot application, we evaluate the returning probability of the three-qutrit system as a function of time:

\begin{align}
p(t) = \Vert\langle 202 | e^{-i t H_{\rm XY}} | 202 \rangle\Vert^2
\label{eq:return_prob}
\end{align}

with the coupling constant $J=0.55$. See the blue curve in Fig. <ref> for the time evolution within $0 \leq t < 5$. Let us consider a Trotterization of (<ref>) with the step size $\theta = 0.025$, resulting in a total of $200$ steps. To employ the above circuit compression strategy, it is crucial to select an appropriate Trotter scheme based on the benchmarking outcomes in Fig. <ref>. We find that for $J=0.55$, choosing $T=T_3$ and $n_b=2$ corresponds to a sufficiently low infidelity,

even less than $10^{-10}.$

We then impose the approximate reflection relation $W_L^{2} \simeq W_R^{2}$ sequentially, by skipping the initial $W_L$ and then substituting the following $W_L^2$ with $W_R^2$. This process repeats by skipping the subsequent $W_L$ and replacing the next $W_L^2$ with $W_R^2$, until no more $W_L^2$ remains. For $200$ Trotter blocks, the substitution is performed $66$ times, leading to the consolidation of $132$ $\mathcal{U}_{x+y}$ gates.
We show the noiseless simulation results from the

approximately compressed circuits as

green triangles in Fig. <ref>.

They demonstrate strong agreement, which is beneficial since it maintains the same level of numerical accuracy while reducing the usage of quantum resources on a qutrit-based quantum computer.

All simulations for exact quantum time dynamics and its Trotterized version shown in Fig. <ref> have been performed using the QuTiP package [31, 32].

§ CONCLUSION AND OUTLOOK

In this paper, we extended the discussion of searching for Yang-Baxter-like turnover relations to qudit-based quantum computing. We explored certain algebraic properties of spin-1 operators and found rigorous Yang-Baxter-like turnover relations for simple qudit models. Regarding more complex qudit models, advanced algebraic relations are challenging to resolve rigorously, but a preliminary and plausible strategy is briefly discussed, based on the conjugation relation between the spin-1 operators. Nevertheless, since large-scale quantum simulation often requires an inexact but sufficiently accurate quantum simulation, we also placed an emphasis on the numerical exploration of advanced relations for the spin-1 system. As a demonstration, we examined the spin-1 XY model and numerically explored advanced circuit turnover relations in the quantum simulation of time dynamics. In particular, we designed a pool of circuit fragment turnover pairs using various Trotterization schemes and numerically examined their fidelity to screen out potential turnover pairs that could be utilized to optimize the deep circuit corresponding to many Trotter steps in the quantum dynamics simulations. Preliminary numerical demonstrations were given on the quantum simulation of the three-qutrit XY model, where the results from our proposed numerical scheme showed great agreement with the exact curve. Remarkably, our numerical scheme can be considered as a prototype of a machine learning process used in circuit compilation and optimization [36] and further development our Yang-Baxter compiler QuYBE [27]. Specifically, the circuit optimization strategy can be boiled down to a combinatorial problem of searching for and performing circuit fragment turnovers in a given circuit with a layered structure, a task that can also be greatly facilitated by high-performance computing (HPC) hardware and brute-force search methods. Further studies in this direction are now underway. Our eventual goal is to develop an efficient parallel Yang-Baxter compiler such as QuYBE [27], which can be used for the compression of a variety of quantum circuits, with the initial target on quantum dynamics circuits.

§ ACKNOWLEDGMENTS

This material is based upon work supported by the  U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Fermilab’s Superconducting Quantum Materials and Systems Center (SQMS) under contract No. DE-AC02-07CH11359 and the National Quantum Center Q-NEXT, which is located in Argonne National Laboratory. O.O. would like to thank SQMS Algorithms Focus for helpful comments.
Y.A. acknowledges support from the U.S. Department of Energy, Office of Science, under contract DE-AC02-06CH11357 at Argonne National Laboratory. We gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. A.B.Ö. thanks Andrey Khesin for insightful discussions on an earlier version of the manuscript. J.K. thanks Matt Reagor for the support on this project and helpful discussions.

§ DERIVATION OF THE TURN-OVER IDENTITIES WITH $\MATHCAL{U}_X$, $\MATHCAL{U}_Y$, AND $\MATHCAL{U}_Z$.

Let $a=x$ and

\begin{align}
\mathcal{U}_x(\alpha) %&= \exp{(i \alpha \tilde{S}^x \otimes \tilde{S}^x)} \notag \\
%&= \mathbf{I}_9 - i\sin(\alpha)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\alpha}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2  \notag \\
&= \begin{pmatrix}
\mathbf{I}_3 & \mathbf{0}_3 & \mathbf{0}_3 \\
\mathbf{0}_3 & \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^x)^2 & \sin(\alpha)\tilde{S}^x \\
\mathbf{0}_3 & -\sin(\alpha)\tilde{S}^x & \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^x)^2
\end{pmatrix} \notag \\ \\
&= \begin{pmatrix}
\mathbf{I}_3 &  \mathbf{0}_3 & \mathbf{0}_3 \\
\mathbf{0}_3 &  \mathbf{A}_{\alpha} & \mathbf{B}_{\alpha} \\
\mathbf{0}_3 & -\mathbf{B}_{\alpha} & \mathbf{A}_{\alpha}
\end{pmatrix}
\end{align}

where we denote

\begin{align}
\mathbf{A}_{\alpha} &= \mathbf{I}_3- 2(\sin(\frac{\alpha}{2})\tilde{S}^x)^2, \\
\mathbf{B}_{\alpha} &= \sin(\alpha)\tilde{S}^x.
\end{align}

Now take the LHS of (<ref>) and replace $\mathcal{R}_n$ with $\mathcal{U}_x$, it is straightforward to show that

\begin{align}
&\big(\mathcal{U}_x(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big)\big(\mathcal{U}_x(\gamma) \otimes \mathbf{I}_3\big) \notag \\
&= \begin{pmatrix}
\mathbf{I}_9 &  \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 &  \mathbf{A}_{\alpha}\otimes\mathbf{I}_3 & \mathbf{B}_{\alpha}\otimes\mathbf{I}_3 \\
\mathbf{0}_9 & -\mathbf{B}_{\alpha}\otimes\mathbf{I}_3 & \mathbf{A}_{\alpha}\otimes\mathbf{I}_3
\end{pmatrix}
\begin{pmatrix}
\mathcal{U}_x(\beta) & \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 & \mathcal{U}_x(\beta) & \mathbf{0}_9\\
\mathbf{0}_9 & \mathbf{0}_9 & \mathcal{U}_x(\beta)
\end{pmatrix}
\begin{pmatrix}
\mathbf{I}_9 &  \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 &  \mathbf{A}_{\gamma}\otimes\mathbf{I}_3 & \mathbf{B}_{\gamma}\otimes\mathbf{I}_3 \\
\mathbf{0}_9 & -\mathbf{B}_{\gamma}\otimes\mathbf{I}_3 & \mathbf{A}_{\gamma}\otimes\mathbf{I}_3
\end{pmatrix} \notag \\
&= \begin{pmatrix}
\mathcal{U}_x(\beta) & \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 & \mathbf{C} & \mathbf{D}\\
\mathbf{0}_9 & -\mathbf{D} & \mathbf{C},
\end{pmatrix}. \label{lhs_x}
\end{align}

\begin{align}
\mathbf{C} &= (\mathbf{A}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{A}_{\gamma}\otimes\mathbf{I}_3) - (\mathbf{B}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{B}_{\gamma}\otimes\mathbf{I}_3), \\
\mathbf{D} &= (\mathbf{A}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{B}_{\gamma}\otimes\mathbf{I}_3) + (\mathbf{B}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{A}_{\gamma}\otimes\mathbf{I}_3)
\end{align}

Similarly, for the RHS of (<ref>), we have

\begin{align}
&\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\delta)\big)\big(\mathcal{U}_x(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\zeta)\big) \notag \\
&= \begin{pmatrix}
\mathcal{U}_x(\delta) & \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 & \mathcal{U}_x(\delta) & \mathbf{0}_9\\
\mathbf{0}_9 & \mathbf{0}_9 & \mathcal{U}_x(\delta)
\end{pmatrix}
\begin{pmatrix}
\mathbf{I}_3 &  \mathbf{0}_3 & \mathbf{0}_3 \\
\mathbf{0}_3 &  \mathbf{A}_{\epsilon} & \mathbf{B}_{\epsilon} \\
\mathbf{0}_3 & -\mathbf{B}_{\epsilon} & \mathbf{A}_{\epsilon}
\end{pmatrix}
\begin{pmatrix}
\mathcal{U}_x(\zeta) & \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 & \mathcal{U}_x(\zeta) & \mathbf{0}_9\\
\mathbf{0}_9 & \mathbf{0}_9 & \mathcal{U}_x(\zeta)
\end{pmatrix} \notag \\
&= \begin{pmatrix}
\mathcal{U}_x(\delta+\zeta) & \mathbf{0}_9 & \mathbf{0}_9 \\
\mathbf{0}_9 & \mathcal{U}_x(\delta)\mathbf{A}_{\epsilon}\mathcal{U}_x(\zeta) & \mathcal{U}_x(\delta)\mathbf{B}_{\epsilon}\mathcal{U}_x(\zeta)\\
\mathbf{0}_9 & -\mathcal{U}_x(\delta)\mathbf{B}_{\epsilon}\mathcal{U}_x(\zeta) & \mathcal{U}_x(\delta)\mathbf{A}_{\epsilon}\mathcal{U}_x(\zeta)
\end{pmatrix}. \label{rhs_x}
\end{align}

Compare (<ref>) with (<ref>), for (<ref>) to hold, the following conditions need to be satisfied,

\begin{align}
&\mathcal{U}_x(\beta) = \mathcal{U}_x(\delta+\zeta); \label{cond_eq1} \\
&\mathbf{C} = \mathcal{U}_x(\delta)\mathbf{A}_{\epsilon}\mathcal{U}_x(\zeta); \label{cond_eq2} \\
&\mathbf{D} = \mathcal{U}_x(\delta)\mathbf{B}_{\epsilon}\mathcal{U}_x(\zeta). \label{cond_eq3}
\end{align}

It's easy to see that (<ref>) is satisfied as long as

\begin{align}
\beta + 2k\pi = \delta+\zeta, ~~ k\in Z. \label{cond1}
\end{align}

In (<ref>)

\begin{align}
&(\mathbf{A}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{A}_{\gamma}\otimes\mathbf{I}_3) \notag \\
&= \Big(\mathbf{I}_9- 2\sin^2(\frac{\alpha}{2}) (\tilde{S}^x)^2\otimes\mathbf{I}_3) \Big) \Big( \mathbf{I}_9 - i\sin(\beta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\beta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \Big(\mathbf{I}_9- 2\sin^2(\frac{\gamma}{2}) (\tilde{S}^x)^2\otimes\mathbf{I}_3) \Big) \notag \\
&= \mathbf{I}_9 + \bigg(\cos(\alpha)\cos(\gamma)-1\bigg)(\tilde{S}^x\otimes\mathbf{I}_3)^2 - i \cos(\alpha)\sin(\beta)\cos(\gamma) (\tilde{S}^x\otimes\tilde{S}^x) - 2 \cos(\alpha)\sin^2(\frac{\beta}{2})\cos(\gamma) (\tilde{S}^x \otimes \tilde{S}^x)^2, \\
&(\mathbf{B}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{B}_{\gamma}\otimes\mathbf{I}_3) \notag \\
&= \Big(\sin(\alpha)(\tilde{S}^x\otimes\mathbf{I}_3) \Big) \Big( \mathbf{I}_9 - i\sin(\beta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\beta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \Big(\sin(\gamma)(\tilde{S}^x\otimes\mathbf{I}_3) \Big) \notag \\
&= \sin(\alpha)\sin(\gamma)(\tilde{S}^x\otimes\mathbf{I}_3)^2 - i \sin(\alpha)\sin(\beta)\sin(\gamma)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin(\alpha)\sin^2(\frac{\beta}{2})\sin(\gamma)(\tilde{S}^x \otimes \tilde{S}^x)^2, \\
\Rightarrow &~~ \mathbf{C} = \mathbf{I}_9 + \bigg(\cos(\alpha+\gamma)-1\bigg)(\tilde{S}^x\otimes\mathbf{I}_3)^2 - i \cos(\alpha+\gamma)\sin(\beta) (\tilde{S}^x \otimes \tilde{S}^x)
- \cos(\alpha+\gamma) \bigg( \cos(\beta) - 1 \bigg) (\tilde{S}^x \otimes \tilde{S}^x)^2, \\
& \mathcal{U}_x(\delta)\mathbf{A}_{\epsilon}\mathcal{U}_x(\zeta) \notag \\
&= \Big( \mathbf{I}_9 - i\sin(\delta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\delta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big)
\Big(\mathbf{I}_9- 2\sin^2(\frac{\epsilon}{2}) (\tilde{S}^x)^2\otimes\mathbf{I}_3) \Big) \notag \\
&~~~~\times \Big( \mathbf{I}_9 - i\sin(\zeta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\zeta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \notag \\
&= \mathbf{I}_9 + \bigg( \cos(\epsilon) - 1 \bigg)(\tilde{S}^x\otimes\mathbf{I}_3)^2 - i \cos(\epsilon)\sin(\delta+\zeta) (\tilde{S}^x \otimes \tilde{S}^x) + \cos(\epsilon)\bigg( \cos(\delta+\zeta) - 1 \bigg)(\tilde{S}^x \otimes \tilde{S}^x)^2.
\end{align}

Similarly, in (<ref>)

\begin{align}
&(\mathbf{A}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{B}_{\gamma}\otimes\mathbf{I}_3) \notag \\
&= \Big(\mathbf{I}_9- 2\sin^2(\frac{\alpha}{2}) (\tilde{S}^x)^2\otimes\mathbf{I}_3) \Big) \Big( \mathbf{I}_9 - i\sin(\beta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\beta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \Big(\sin(\gamma)(\tilde{S}^x\otimes\mathbf{I}_3) \Big) \notag \\
&= \cos(\alpha)\sin(\gamma)(\tilde{S}^x\otimes\mathbf{I}_3) -i \cos(\alpha)\sin(\beta)\sin(\gamma)\bigg( (\tilde{S}^x)^2 \otimes \tilde{S}^x\bigg) - 2\cos(\alpha)\sin^2(\frac{\beta}{2})\sin(\gamma)\bigg( \tilde{S}^x \otimes (\tilde{S}^x)^2 \bigg)\\
&(\mathbf{B}_{\alpha}\otimes\mathbf{I}_3)\mathcal{U}_x(\beta)(\mathbf{A}_{\gamma}\otimes\mathbf{I}_3) \notag \\
&= \Big(\sin(\alpha)(\tilde{S}^x\otimes\mathbf{I}_3) \Big) \Big( \mathbf{I}_9 - i\sin(\beta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\beta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \Big(\mathbf{I}_9- 2\sin^2(\frac{\gamma}{2}) (\tilde{S}^x)^2\otimes\mathbf{I}_3) \Big) \notag \\
&= \sin(\alpha)\cos(\gamma)(\tilde{S}^x\otimes\mathbf{I}_3) -i \sin(\alpha)\sin(\beta)\cos(\gamma)\bigg( (\tilde{S}^x)^2 \otimes \tilde{S}^x\bigg) - 2\sin(\alpha)\sin^2(\frac{\beta}{2})\cos(\gamma)\bigg( \tilde{S}^x \otimes (\tilde{S}^x)^2 \bigg)\\
\Rightarrow &~~ \mathbf{D} = \sin(\alpha+\gamma)(\tilde{S}^x\otimes\mathbf{I}_3) -i \sin(\alpha+\gamma)\sin(\beta)\bigg( (\tilde{S}^x)^2 \otimes \tilde{S}^x\bigg) + \sin(\alpha+\gamma)\bigg(\cos(\beta) - 1 \bigg)\bigg( \tilde{S}^x \otimes (\tilde{S}^x)^2 \bigg)\\
& \mathcal{U}_x(\delta)\mathbf{B}_{\epsilon}\mathcal{U}_x(\zeta) \notag \\
&= \Big( \mathbf{I}_9 - i\sin(\delta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\delta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big)     \Big(\sin(\epsilon) (\tilde{S}^x\otimes\mathbf{I}_3) \Big) \Big( \mathbf{I}_9 - i\sin(\zeta)(\tilde{S}^x \otimes \tilde{S}^x) - 2\sin^2(\frac{\zeta}{2})(\tilde{S}^x \otimes \tilde{S}^x)^2 \Big) \notag \\
&= \sin(\epsilon)(\tilde{S}^x\otimes\mathbf{I}_3) -i \sin(\epsilon)\sin(\delta+\zeta)\bigg( (\tilde{S}^x)^2 \otimes \tilde{S}^x\bigg) + \sin(\epsilon)\bigg(\cos(\delta+\zeta)-1\bigg)\bigg( \tilde{S}^x \otimes (\tilde{S}^x)^2 \bigg).
\end{align}

After a term-by-term comparison, one can see that the following conditions need to be satisfied

\begin{align}
\bigg\{\begin{array}{l}
\alpha + \gamma = \epsilon + 2k\pi  \\
\delta + \zeta = \beta + 2k\pi
\end{array},
~~ k\in Z.
\end{align}

for the following equation to hold

\begin{align}
&\big(\mathcal{U}_x(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big)\big(\mathcal{U}_x(\gamma) \otimes \mathbf{I}_3\big) \notag \\
&~~~~~~~~=\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\delta)\big)\big(\mathcal{U}_x(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\zeta)\big).
\end{align}

From (<ref>) a further proof can be obtained for $a=y$

\begin{align}
&\big(\mathcal{U}_y(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_y(\beta)\big)\big(\mathcal{U}_y(\gamma) \otimes \mathbf{I}_3\big) \notag \\
&= \bigg( (P_y\otimes P_y \otimes \mathbf{I}_3) \big(\mathcal{U}_x(\alpha) \otimes \mathbf{I}_3\big) (P_y\otimes P_y \otimes \mathbf{I}_3) \bigg)
\bigg( (\mathbf{I}_3 \otimes P_y\otimes P_y ) \big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big) (\mathbf{I}_3 \otimes P_y\otimes P_y ) \bigg) \notag \\
&~~~~ \times \bigg( (P_y\otimes P_y \otimes \mathbf{I}_3) \big(\mathcal{U}_x(\gamma) \otimes \mathbf{I}_3\big) (P_y\otimes P_y \otimes \mathbf{I}_3) \bigg) \notag \\
&= (P_y\otimes P_y \otimes P_y) \big(\mathcal{U}_x(\alpha) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big)\big(\mathcal{U}_x(\gamma) \otimes \mathbf{I}_3\big) (P_y\otimes P_y \otimes P_y) \notag \\
&= (P_y\otimes P_y \otimes P_y) \big(\mathbf{I}_3 \otimes \mathcal{U}_x(\delta)\big)\big(\mathcal{U}_x(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_x(\zeta)\big) (P_y\otimes P_y \otimes P_y) \notag \\
&= \bigg( (\mathbf{I}_3 \otimes P_y\otimes P_y ) \big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big) (\mathbf{I}_3 \otimes P_y\otimes P_y ) \bigg) \bigg( (P_y\otimes P_y \otimes \mathbf{I}_3) \big(\mathcal{U}_x(\alpha) \otimes \mathbf{I}_3\big) (P_y\otimes P_y \otimes \mathbf{I}_3) \bigg) \notag \\
&~~~~ \bigg( (\mathbf{I}_3 \otimes P_y\otimes P_y ) \big(\mathbf{I}_3 \otimes \mathcal{U}_x(\beta)\big) (\mathbf{I}_3 \otimes P_y\otimes P_y ) \bigg) \notag \\
&= \big(\mathbf{I}_3 \otimes \mathcal{U}_y(\delta)\big)\big(\mathcal{U}_y(\epsilon) \otimes \mathbf{I}_3\big)\big(\mathbf{I}_3 \otimes \mathcal{U}_y(\zeta)\big).
\end{align}

Similar process can be followed from (<ref>) for the proof for $a=z$.

§ QUTRIT IDENTITIES FROM QUBIT IDENTITIES

The adjoint spin-1 matrices $\tilde{S}^a$ can be seen as a $2\times 2$ block embedding of Pauli-$Y$ into a $3\times 3$ matrix, implying that the unitaries $\mathcal{U}_a$ primarily affect only two levels of the qutrit. We can build a permutation $P$ that separates the spectator levels from the levels actively involved in the unitary operations.

Taking $\tilde{S}_x$ as an example, from (<ref>) we have

\begin{align*}
P =
\resizebox{0.9\hsize}{!}{$\begin{psmallmatrix}
4 & 5 & 6  & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23 & 24 \\
5 & 6 & 4  & 11 & 12 & 15 & 19 & 20 & 16 & 21 & 22 & 10 & 13 & 14 & 17 & 23 & 24 & 18
\end{psmallmatrix}$}
\end{align*}

which rearranges all 3-qutrit states into the following direct sum between subspaces: $|000\rangle  \oplus |00a\rangle  \oplus |0a0\rangle \oplus |0ab\rangle \oplus |a00\rangle  \oplus |a0b\rangle \oplus |ab0\rangle \oplus |abc\rangle $ with $a,b,c \in \{1,2\}$. Such permutation decomposes the 2-qutrit unitary action $\mathcal{U}_x$ into the direct sum $\mathbf{1}_{1} \oplus \exp(-i \alpha Y_1) \oplus \exp(-i \alpha Y_2)\oplus \exp(-i \alpha Y_1 \otimes Y_2)$ where $Y_{1,2}$ act on effective “qubits" obtained from restricting qutrits onto the two-levels $|1\rangle$, $|2\rangle$.

Applying the above permutation to the qutrit turnover relation for $\mathcal{U}_x$,
both sides of (<ref>) takes a block-diagonal form corresponding to some circuit identities on effective two-level systems.

* $|000\rangle$: The identity $1=1$ holds trivially.

* $|00a\rangle$: $e^{-i \alpha Y_3} = e^{-i \delta Y_3}  e^{-i \zeta Y_3} $

* $|0a0\rangle$: $e^{-i \alpha Y_2}e^{-i \beta Y_2}e^{-i \gamma Y_2} = e^{-i\delta Y_2} e^{-i \epsilon Y_2} e^{-i \zeta Y_2}$

* $|a00\rangle$: $e^{-i \alpha Y_1}e^{-i \gamma Y_1} = e^{-i \epsilon Y_1}$

* $|a0b\rangle$: $e^{-i \alpha Y_1}e^{-i \beta Y_3}e^{-i \gamma Y_1} = e^{-i \delta Y_3}e^{-i \epsilon Y_1}e^{-i \zeta Y_3}$

* $|0ab\rangle$

\begin{align*}
e^{-i \alpha Y_2}e^{-i \beta Y_2 \otimes Y_3}e^{-i \gamma Y_2} = e^{-i\delta Y_2 \otimes Y_3} e^{-i \epsilon Y_2} e^{-i \zeta Y_2 \otimes Y_3}
\end{align*}

* $|ab0\rangle$

\begin{align*}
e^{-i\alpha Y_1 \otimes Y_2} e^{-i \beta Y_2} e^{-i \gamma Y_1 \otimes Y_2} = e^{-i \delta Y_2}e^{-i \epsilon Y_1 \otimes Y_2}e^{-i \zeta Y_2}
\end{align*}

* $|abc\rangle$: It returns the qubit turnover relation of [46].

The relations on $|00a\rangle$, $|0a0\rangle$, $|a00\rangle$, $|a0b\rangle$ subspaces are obviously true if the circuit parameters satisfy  (<ref>). The $|0ab\rangle$ and $|ab0\rangle$ relations are also straightforward to verify since $[e^{-i  u Y_i}, \,e^{-i v Y_i \otimes Y_j}] = 0$. The only remaining non-trivial circuit relation comes from the subspace $|abc\rangle$, equivalent to the “qubit" turnover relation shown in [46].

[1]
Yosi Atia and Dorit Aharonov.
Fast-forwarding of hamiltonians and exponentially precise
Nat. Commun., 8(1):1–9, 2017.

[2]
Fergus Barratt, James Dborin, Matthias Bal, Vid Stojevic, Frank Pollmann, and
Andrew G Green.
Parallel quantum simulation of large systems on small nisq computers.
npj Quantum Inf., 7(1):1–7, 2021.

[3]
Murray T Batchelor and Angela Foerster.
Yang–Baxter integrable models in experiments: from condensed
matter to ultracold atoms.
J. Phys. A: Math. Theor., 49(17):173001, 2016.

[4]
Rodney J Baxter.
Exactly solved models in statistical mechanics.
Elsevier, 2016.

[5]
Helle Bechmann-Pasquinucci and Asher Peres.
Quantum cryptography with 3-state systems.
Phys. Rev. Lett., 85:3313–3316, Oct 2000.

[6]
Noah F Berthusen, Thaís V Trevisan, Thomas Iadecola, and Peter P Orth.
Quantum dynamics simulations beyond the coherence time on noisy
intermediate-scale quantum hardware by variational trotter compression.
Physical Review Research, 4(2):023097, 2022.

[7]
M. S. Blok, V. V. Ramasesh, T. Schuster, K. O'Brien, J. M. Kreikebaum,
D. Dahlen, A. Morvan, B. Yoshida, N. Y. Yao, and I. Siddiqi.
Quantum information scrambling on a superconducting qutrit processor.
Phys. Rev. X, 11:021010, Apr 2021.

[8]
Alex Bocharov, Martin Roetteler, and Krysta M. Svore.
Factoring with qutrits: Shor's algorithm on ternary and metaplectic
quantum architectures.
Phys. Rev. A, 96:012306, Jul 2017.

[9]
S. B. Bravyi and A. Yu. Kitaev.
Quantum codes on a lattice with boundary.
ArXiv Preprint, 1998.

[10]
Gavin K. Brennen, Dianne P. O'Leary, and Stephen S. Bullock.
Criteria for exact qudit universality.
Phys. Rev. A, 71:052318, May 2005.

[11]
D. Bruß and C. Macchiavello.
Optimal eavesdropping in cryptography with three-dimensional quantum
Phys. Rev. Lett., 88:127901, Mar 2002.

[12]
Stephen S. Bullock, Dianne P. O'Leary, and Gavin K. Brennen.
Asymptotically optimal quantum circuits for $d$-level systems.
Phys. Rev. Lett., 94:230502, Jun 2005.

[13]
Earl T. Campbell.
Enhanced fault-tolerant quantum computing in $d$-level systems.
Phys. Rev. Lett., 113:230501, Dec 2014.

[14]
Earl T. Campbell, Hussain Anwar, and Dan E. Browne.
Magic-state distillation in all prime dimensions using quantum
reed-muller codes.
Phys. Rev. X, 2:041021, Dec 2012.

[15]
Daan Camps, Efekan Kökcü, Lindsay Bassman Oftelie, Wibe A. de Jong,
Alexander F. Kemper, and Roel Van Beeumen.
An algebraic quantum circuit compression algorithm for hamiltonian
SIAM J. Matrix Anal. Appl., 43(3):1084–1108, 2022.

[16]
Rebecca S. Chen.
Generalized yang–baxter equations and braiding quantum gates.
J. Knot Theory Ramif., 21(09):1250087, 2012.

[17]
Zijun Chen, Kevin J. Satzinger, Juan Atalaya, Alexander N. Korotkov, Andrew
Dunsworth, Daniel Sank, Chris Quintana, Matt McEwen, Rami Barends, Paul V.
Klimov, and et al.
Exponential suppression of bit or phase errors with cyclic error
Nature, 595(7867):383–387, jul 2021.

[18]
Yulin Chi, Jieshan Huang, Zhanchuan Zhang, Jun Mao, Zinan Zhou, Xiaojiong Chen,
Chonghao Zhai, Jueming Bao, Tianxiang Dai, Huihong Yuan, and et al.
A programmable qudit-based quantum processor.
Nat. Commun., 13:1166, mar 2022.

[19]
Cristina Cirstoiu, Zoe Holmes, Joseph Iosue, Lukasz Cincio, Patrick J Coles,
and Andrew Sornborger.
Variational fast forwarding for quantum simulation beyond the
coherence time.
npj Quantum Information, 6(1):1–10, 2020.

[20]
Eric Dennis, Alexei Kitaev, Andrew Landahl, and John Preskill.
Topological quantum memory.
J. Math. Phys., 43(9):4452–4505, 08 2002.

[21]
David P. DiVincenzo.
Two-bit gates are universal for quantum computation.
Phys. Rev. A, 51:1015–1022, Feb 1995.

[22]
Laurin E. Fischer, Alessandro Chiesa, Francesco Tacchino, Daniel J. Egger,
Stefano Carretta, and Ivano Tavernelli.
Towards universal gate synthesis and error correction in transmon
ArXiv Preprint, page arXiv:2212.04496, 2022.

[23]
César Galindo, Seung-Moon Hong, and Eric C. Rowell.
Generalized and quasi-localizations of braid group representations.
ArXiv Preprint, page arXiv:1105.5048, 2011.

[24]
Mo-Lin Ge, Kang Xue, Ruo-Yang Zhang, and Qing Zhao.
Yang–Baxter equations and quantum entanglements.
Quantum Inf. Process., 15(12):5211–5242, 2016.

[25]
Pranav Gokhale, Jonathan M. Baker, Casey Duckering, Natalie C. Brown,
Kenneth R. Brown, and Frederic T. Chong.
Asymptotic improvements to quantum circuits via qutrits.
In Proceedings of the 46th International Symposium on Computer
Architecture, ISCA '19, page 554–566, New York, NY, USA, 2019. Association
for Computing Machinery.

[26]
Noah Goss, Alexis Morvan, Brian Marinelli, Bradley K. Mitchell, Long B. Nguyen,
Ravi K. Naik, Larry Chen, Christian Jünger, John Mark Kreikebaum, David I.
Santiago, Joel J. Wallman, and Irfan Siddiqi.
High-fidelity qutrit entangling gates for superconducting circuits.
Nat. Commun., 13(1), dec 2022.

[27]
Sahil Gulania, Zichang He, Bo Peng, Niranjan Govind, and Yuri Alexeev.
QuYBE - an algebraic compiler for quantum circuit compression.
In 2022 IEEE/ACM 7th Symposium on Edge Computing (SEC), pages
406–410. IEEE, 2022.

[28]
Erik Gustafson.
Noise improvements in quantum simulations of sqed using qutrits.
ArXiv Preprint, page arXiv:2201.04546, 2022.

[29]
Pavel Hrmo, Benjamin Wilhelm, Lukas Gerster, Martin W. van Mourik, Marcus
Huber, Rainer Blatt, Philipp Schindler, Thomas Monz, and Martin Ringbauer.
Native qudit entanglement in a trapped ion quantum processor.
ArXiv Preprint, page arXiv:2206.04104, 2022.

[30]
A. P. Isaev.
Lectures on quantum groups and yang-baxter equations.
ArXiv Preprint, page arXiv:2206.08902, 2022.

[31]
J Robert Johansson, Paul D Nation, and Franco Nori.
Qutip: An open-source python framework for the dynamics of open
quantum systems.
Computer Physics Communications, 183(8):1760–1772, 2012.

[32]
J.R. Johansson, P.D. Nation, and Franco Nori.
Qutip 2: A python framework for the dynamics of open quantum systems.
Computer Physics Communications, 184(4):1234–1240, 2013.

[33]
Eliot Kapit.
Hardware-efficient and fully autonomous quantum error correction in
superconducting circuits.
Phys. Rev. Lett., 116:150501, Apr 2016.

[34]
L. H. Kauffman and S. J. Jr. Lomonaco.
Topological quantum information theory.
In S. J. Lomonaco, editor, Proceedings of Symposia in Applied
Mathematics, Vol. 68, Washington, DC, 2010. AMS.

[35]
Efekan Kökcü, Daan Camps, Lindsay Bassman, J. K. Freericks, Wibe A.
de Jong, Roel Van Beeumen, and Alexander F. Kemper.
Algebraic compression of quantum circuits for Hamiltonian
Phys. Rev. A, 105(3), mar 2022.

[36]
Mario Krenn, Jonas Landgraf, Thomas Foesel, and Florian Marquardt.
Artificial intelligence and machine learning for quantum
Phys. Rev. A, 107:010101, Jan 2023.

[37]
Sebastian Krinner, Nathan Lacroix, Ants Remm, Agustin Di Paolo, Elie Genois,
Catherine Leroux, Christoph Hellings, Stefania Lazar, Francois Swiadek,
Johannes Herrmann, and et al.
Realizing repeated quantum error correction in a distance-three
surface code.
Nature, 605(7911):669–674, may 2022.

[38]
Sheng-Hsuan Lin, Rohit Dilip, Andrew G Green, Adam Smith, and Frank Pollmann.
Real-and imaginary-time evolution with compressed quantum circuits.
PRX Quantum, 2(1):010342, 2021.

[39]
Ritajit Majumdar, Saikat Basu, Shibashis Ghosh, and Susmita Sur-Kolay.
Quantum error-correcting code for ternary logic.
Phys. Rev. A, 97:052302, May 2018.

[40]
J. F. Marques, B. M. Varbanov, M. S. Moreira, H. Ali, N. Muthusubramanian,
C. Zachariadis, F. Battistel, M. Beekman, N. Haider, W. Vlothuizen, and
et al.
Logical-qubit operations in an error-detecting surface code.
Nat. Phys., 18(1):80–86, dec 2021.

[41]
Sreraman Muralidharan, Chang-Ling Zou, Linshu Li, Jianming Wen, and Liang
Overcoming erasure errors with multilevel systems.
New J. Phys., 19(1):013026, jan 2017.

[42]
Ashok Muthukrishnan and C. R. Stroud.
Multivalued logic gates for quantum computation.
Phys. Rev. A, 62:052309, Oct 2000.

[43]
Chetan Nayak, Steven H. Simon, Ady Stern, Michael Freedman, and Sankar
Das Sarma.
Non-Abelian anyons and topological quantum computation.
Rev. Mod. Phys., 80:1083–1159, Sep 2008.

[44]
Lindsay Bassman Oftelie, Roel Van Beeumen, Ed Younis, Ethan Smith, Costin
Iancu, and Wibe A. de Jong.
Constant-depth circuits for dynamic simulations of materials on
quantum computers.
Mater. Theory, 6, 2022.

[45]
Archimedes Pavlidis and Emmanuel Floratos.
Quantum-fourier-transform-based quantum arithmetic with qudits.
Phys. Rev. A, 103:032417, Mar 2021.

[46]
Bo Peng, Sahil Gulania, Yuri Alexeev, and Niranjan Govind.
Quantum time dynamics employing the yang-baxter equation for circuit
Phys. Rev. A, 106:012412, Jul 2022.

[47]
Martin Ringbauer, Michael Meth, Lukas Postler, Roman Stricker, Rainer Blatt,
Philipp Schindler, and Thomas Monz.
A universal qudit quantum processor with trapped ions.
Nat. Phys., 18(9):1053–1057, jul 2022.

[48]
Eric C. Rowell.
A quaternionic braid representation (after goldschmidt and jones).
ArXiv Preprint, page arXiv:1006.4808, 2010.

[49]
Eric C. Rowell, Yong Zhang, Yong-Shi Wu, and Mo-Lin Ge.
Extraspecial two-groups, generalized yang-baxter equations and
braiding quantum gates.
Quantum Inf. Comput., 10(7&8):685–702, 2010.

[50]
Alipasha Vaziri, Gregor Weihs, and Anton Zeilinger.
Experimental two-photon, three-dimensional entanglement for quantum
Phys. Rev. Lett., 89:240401, Nov 2002.

[51]
F Anvari Vind, Angela Foerster, Ivan S Oliveira, Roberto Silva Sarthour, Diogo
de Oliveira Soares-Pinto, Alexandre Martins de Souza, and Itzhak Roditi.
Experimental realization of the Yang–Baxter equation via NMR
Scientific reports, 6(1):1–8, 2016.

[52]
Chen-Ning Yang.
Some exact results for the many-body problem in one dimension with
repulsive delta-function interaction.
Phys. Rev. Lett., 19(23):1312, 1967.

[53]
Y. Zhang.
Integrable quantum computation.
Quantum Inf. Process., 12:631–639, 2013.

[54]
Youwei Zhao, Yangsen Ye, He-Liang Huang, Yiming Zhang, Dachao Wu, Huijie Guan,
Qingling Zhu, Zuolin Wei, Tan He, Sirui Cao, and et al.
Realization of an error-correcting surface code with superconducting
Phys. Rev. Lett., 129:030501, Jul 2022.

[55]
D. L. Zhou, B. Zeng, Z. Xu, and C. P. Sun.
Quantum computation based on d-level cluster state.
Phys. Rev. A, 68:062303, Dec 2003.