Quantum time dynamics of 1D-Heisenberg models employing the Yang-Baxter equation for circuit compression

Sahil Gulania sgulania@anl.gov Mathematics and Computer Science Division, Argonne National Laboratory, Lemont, Illinois 60439, United States Bo Peng peng398@pnnl.gov Physical and Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, Washington 99352, United States Yuri Alexeev yuri@anl.gov Computational Science Division, Argonne National Laboratory, Lemont, Illinois 60439, United States Niranjan Govind niri.govind@pnnl.gov Physical and Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, Washington 99352, United States[]

Abstract

Quantum time dynamics (QTD) is considered a promising problem for quantum supremacy on near-term quantum computers. However, QTD quantum circuits grow with increasing time simulations. This study focuses on simulating the time dynamics of 1-D integrable spin chains with nearest neighbor interactions. We show how the quantum Yang-Baxter equation can be exploited to compress and produce a shallow quantum circuit. With this compression scheme, the depth of the quantum circuit becomes independent of step size and only depends on the number of spins. We show that the compressed circuit scales quadratically with system size, which allows for the simulations of time dynamics of very large 1-D spin chains. We derive the compressed circuit representations for different special cases of the Heisenberg Hamiltonian. We compare and demonstrate the effectiveness of this approach by performing simulations on quantum computers.

^†^†preprint: AIP/123-QED

I Introduction

Simulation of statistical mechanical models is a vital application for classical and quantum computing [1]. It is also well known that the partition function of a $d+1$ -dimensional classical system can be mapped to the partition function of a $d$ -dimensional quantum system [2, 3, 4]. This deep classical to quantum connection can give insights into quantum universality classes using appropriate classical counterparts. It might allow one to use these effective models to study critical points and phase transitions of magnetic systems, where the spins of the magnetic systems are treated quantum mechanically.

Quantum Ising and Heisenberg models [5] represent some of the simplest models that can describe the behavior of magnetic systems. However, 2D and 3D-quantum lattice simulations remain challenging for classical computing. A possible solution is to use quantum computing, given that it is only natural to simulate quantum systems with quantum computers as suggested by Benioff [6] and Feynman [7]. The idea of using quantum computing on quantum devices is compelling since it paves the way for systematic improvements of quantum technologies. The exploration of this approach is at the core of several research efforts in quantum information sciences [8]. In particular, the development of new generations of quantum models and associated simulations on existing and upcoming Noisy Intermediate-Scale Quantum (NISQ) [9] devices is of particular interest.

Quantum integrable systems [10] typically refer to systems where the dynamics are two-body reducible. Put another way, even though the Hilbert space increases exponentially with increasing system size in these systems, the two-body reducibility, in combination with the algebraic Bethe ansatz, [11, 12, 5] can be used to obtain explicit solutions, under certain conditions, by solving a set of non-linear equations that scales only linearly with system size. The quantum Yang–Baxter equation or star-triangle relation [13, 14, 10, 15] is a consequence of this factorization. The algebraic formulation of quantum integrable systems makes them ideal tools to study a broad range of low dimension physical models. Historically, the isotropic interacting quantum spin chain, or Heisenberg model [16] was the first quantum integrable system, whose exact eigenstates were obtained through the Bethe ansatz approach as a superposition of plane-waves [11, 17, 18, 19, 20, 21]. Other quantum integrable models include the Lieb–Liniger model [22, 23], the Hubbard model [24], Calogero-Sutherland model [25, 26, 27, 28, 29, 5], models from quantum field theory like the sine-Gordon model [30, 31], and several sub-classes of the Heisenberg model (e.g. XXZ model) [16, 32].

An important yet open question in quantum integrable systems is time evolution, or the response dynamics where the system responds to a change of parameters. This requires sufficiently accurate control of the time evolution of the system which is governed by the time-dependent Schrödinger or Dirac equation of the quantum state in the Hilbert space. The time evolution problem is also closely related to the computation of the asymptotic state of quantum integrable models, or in the more general sense the thermalization and ergodic/non-ergodic behaviors of these models. Several conjectures have been proposed in this context [33, 34, 35]. For example, in the Heisenberg spin chain, quantum quenches of the XXZ model has been studied by embedding the generalized Gibbs ensemble hypothesis into the quantum transfer matrix framework [36, 37]. Nevertheless, these studies are still far from conclusive. Other types of quantum integrability are known in explicitly time-dependent quantum problems, such as the driven Tavis-Cummings model [38, 39].

It is well known that quantum circuits representing quantum time dynamics (QTD) grow with increasing time simulations. In the present era of noisy quantum computers, it has become a necessity for circuits to be as shallow as possible for meaningful results. With this as the overarching theme, we focus on QTD of 1-D integrable spin chains with nearest neighbor interactions in this paper. We show how the quantum Yang-Baxter equation can be used to compress and produce a shallow quantum circuit, where the depth becomes independent of step size and only depends on the number of spins. The compressed circuit scales quadratically with system size, which allows for the simulations of time dynamics of very large 1-D spin chains. Compressed circuit representations are derived for different special cases of the Heisenberg Hamiltonian. As a proof of principle, we demonstrate the effectiveness of this approach by performing simulations of the Heisenberg XY model (or XY model in brief) on quantum devices. The time evolution of the XY model is an active area of research and has been approached from many directions. Verstraete and co-workers [40] utilized the quantum Fourier transform with the Bogoliubov transformation to perform time dynamics efficiently on a quantum computer. Recently, in a similar study to this work, Bassman and co-workers [41] also reported simulations of the XY model by conjecturing the relationship between the reflection symmetry and YBE and provided numerical evidence for the YBE transformation. In this paper we go further by putting the relationship on a firm footing and derive the analytical expressions. For the rest of this paper we will tacitly assume the quantum Yang-Baxter equation (YBE) and omit quantum for brevity.

II Theory

For completeness, we start by reviewing relevant background material, with a brief introduction to the Heisenberg Hamiltonian, quantum time dynamics, and the Yang-Baxter equation.

II.1 Heisenberg Hamiltonian

The Heisenberg Hamiltonian [42, 43, 44] is widely used to study magnetic systems, where the magnetic spins are treated quantum mechanically. The Hamiltonian, including only spin-spin interactions, can be written as

\hat{H}=-\sum_{\alpha}\{J_{\alpha}\sum_{i=1}^{N-1}\sigma_{i}^{\alpha}\otimes\sigma_{i+1}^{\alpha}\}

(1)

where, $\alpha$ sums over $\{x,y,z\}$ , the coupling parameters $J_{\alpha}$ denotes the exchange interaction between nearest neighbour spins along the $\alpha-$ direction, $\sigma^{\alpha}_{i}$ is the $\alpha$ -Pauli operator on the $i$ -th spin. Interaction with the magnetic field can be included in the above Hamiltonian as

\hat{H}_{in}(t)=\hat{H}-h_{\beta}(t)\sum_{i=1}^{N}\sigma_{i}^{\beta}

(2)

where, $h_{\beta}(t)$ is the time amplitude of the external magnetic field along the $\beta\in\{x,y,z\}$ direction. Several variations of this model are known in the literature, which are categorized depending on the relation between $J_{x}$ , $J_{y}$ , and $J_{z}$ . This Heisenberg Hamiltonian represents quantities based on the electronic structure of the system, where the Coulomb interaction and hopping are mapped onto spin variables [45, 46].

A simple variant of the Heisenberg model is the one-dimensional XY model that was first introduced and solved by Lieb, Schultz, and Mattis [47] in the absence of a magnetic field, and later by Katsura [48, 49] and Niemeijer [50] in a finite external field. The XY model describes a one-dimensional lattice with spin variables labeling every lattice point. The spins are limited to interact only with their nearest neighbors in an anisotropic way. As a quantum integrable model, all static correlation functions of the XY model can be computed in terms of Toeplitz determinants [51]. For example, beside the fundamental correlation functions, one can also compute the Emptiness Formation Probability [52, 53, 54], the Von Neumann and Renyi entanglement entropies [55, 56, 57, 58, 59, 60], and even some non-equilibrium properties [61, 62, 63, 64, 65, 66, 67] of the XY model.

It has been shown [68, 69, 70, 71] that the excitations in the XY model are non-local free fermions that give rise to a non-trivial phase diagram at zero temperature. There exists two quantum phase transitions (QPTs) in this diagram featuring the universality of the anti-ferromagnetic Heisenberg chain (i.e. the isotropic XX model) and the quantum Ising model, respectively. In particular, the isotropic XX model corresponds to free fermions hopping on a lattice, and the Ising model corresponds to a transition from a doubly degenerate ground state to a single ground state [72, 73].

II.2 Time Evolution

Quantum state evolution [74, 75] is governed by the Schr $\ddot{o}$ dinger or Dirac equation

i\hbar\frac{\partial}{\partial t}\ket{\psi(t)}=\hat{H}\ket{\psi(t)}

(3)

The solution to the above equation can be expressed as

\displaystyle\ket{\psi(t)}=e^{-i\hat{H}t/\hbar}\ket{\psi(0)}

(4)

where, $e^{-i\hat{H}t/\hbar}$ is the evolution operator. In the 1D-Heisenberg model, with the exception of $N=2$ , all the elements in the Hamiltonian do not commute with each other, and hence the exponential of $\hat{H}$ cannot be written as a product of exponentials. For $N=2$ ,

e^{-i\hat{H}t/\hbar}=\prod_{\alpha}e^{iJ_{\alpha}t(\sigma_{1}^{\alpha}\otimes\sigma_{2}^{\alpha})/\hbar}

(5)

where each term is straightforward to evaluate as shown below

$\displaystyle e^{iJ_{x}t(\sigma_{1}^{x}\otimes\sigma_{2}^{x})/\hbar}=$	$\displaystyle\begin{pmatrix}\cos(\theta_{x})&0&0&i\sin(\theta_{x})\\ 0&\cos(\theta_{x})&i\sin(\theta_{x})&0\\ 0&i\sin(\theta_{x})&\cos(\theta_{x})&0\\ i\sin(\theta_{x})&0&0&\cos(\theta_{x})\end{pmatrix}$	(6)

$\displaystyle e^{iJ_{y}t(\sigma_{1}^{y}\otimes\sigma_{2}^{y})/\hbar}=$	$\displaystyle\begin{pmatrix}\cos(\theta_{y})&0&0&-i\sin(\theta_{y})\\ 0&\cos(\theta_{y})&i\sin(\theta_{y})&0\\ 0&i\sin(\theta_{y})&\cos(\theta_{y})&0\\ -i\sin(\theta_{y})&0&0&\cos(\theta_{y})\end{pmatrix}$

$\displaystyle e^{iJ_{z}t(\sigma_{1}^{z}\otimes\sigma_{2}^{z})/\hbar}=$	$\displaystyle\begin{pmatrix}e^{i\theta_{z}}&0&0&0\\ 0&e^{-i\theta_{z}}&0&0\\ 0&0&e^{-i\theta_{z}}&0\\ 0&0&0&e^{i\theta_{z}}\end{pmatrix}$

where, $\theta_{\alpha}=tJ_{\alpha}/\hbar$ . For $N=3$ , there are terms that do not commute. For instance, if $p_{12}$ represents the Heisenberg interaction (Eq. (5)) between spins 1 and 2 and $p_{23}$ represents interaction between spins 2 and 3, then $p_{12}$ does not commute with $p_{23}$ . As a result, one cannot decompose the time evolution operator as a product of two-body evolution operators. The Trotter decomposition can be used to rewrite the time evolution operator in terms of two-body components as follows,

	$\displaystyle e^{-i\hat{H}t/\hbar}=$	$\displaystyle\big{[}\big{(}\prod_{\alpha}e^{i\theta_{\alpha}(\sigma_{1}^{\alpha}\otimes\sigma_{2}^{\alpha}\otimes\mathbb{1})/n}\big{)}\times\big{(}\prod_{\alpha}e^{i\theta_{\alpha}(\mathbb{1}\otimes\sigma_{2}^{\alpha}\otimes\sigma_{3}^{\alpha})/n}\big{)}\big{]}^{n}$		(7)
		$\displaystyle+\mathcal{O}(t/n)$		(7)

where the error scales linearly with time step i.e. $t/n$ , which can be a significant source of error. This can be mitigated by taking a smaller step size. However, this results in an overall increase in the computation cost. Therefore, there is a need to balance accuracy and computation cost.

Extending the evolution operator to systems with $N>3$ , one can observed that there are two major commuting families as shown in Fig.1. All elements in the orange family commute and all elements in blue family commute. Therefore the evolution operator can be written as a product of exponentials within the families without Trotter decomposition.

Refer to caption — Figure 1: 1D-spin chain showing two families (orange and blue) spin interaction which commute within family.

II.3 Yang-Baxter Equation

The Yang-Baxter equation (YBE) was introduced independently in theoretical physics by Yang [13] in the late 1960s, and by Baxter [76] in statistical mechanics in the early 1970s. This relation has also received much attention in many areas of theoretical physics, classification of knots, scattering of subatomic particles, nuclear magnetic resonance, ultracold atoms, and more recently in quantum information science [77, 78, 79, 80, 81, 82].

The YBE connection to quantum computing originates from the pursuit of figuring out the relationship between topological entanglement, quantum entanglement, and quantum computational universality. In particular, how the global topological relationship in spaces (e.g. knotting and linking) corresponds to the entangled quantum states, and how the CNOT gate, for instance, can in turn be replaced by another unitary gate $R$ to maintain universality. It turns out these unitary $R$ gates that serve to maintain the universality of quantum computation and also as solutions for the condition of topological braiding are unitary solutions to the YBEs [14] Briefly, the relation is a consistency or exchange condition that allows one to factorize the interactions of three bodies into a sequence of pairwise interactions under certain conditions. Formally, this can be written as

(\mathcal{R}\otimes\mathbb{1})(\mathbb{1}\otimes\mathcal{R})(\mathcal{R}\otimes\mathbb{1})=(\mathbb{1}\otimes\mathcal{R})(\mathcal{R}\otimes\mathbb{1})(\mathbb{1}\otimes\mathcal{R})

(8)

where the $\mathcal{R}$ operator is a linear mapping $\mathcal{R}:V\otimes V\rightarrow V\otimes V$ defined as a two-fold tensor product generalizing the permutation of vector space $V$ . This relation also yields a sufficiency condition for quantum integrability in one dimensional quantum systems, and provides a systematic approach to construct integrable models. Since a detailed discussion of this topic is beyond the scope of this paper, we refer the reader to more comprehensive works and reviews on the subject [10, 14].

III Circuit Representation of the Time Evolution Operator

Since the evolution operator is a unitary matrix, there exists a quantum circuit that can perform this operation efficiently on a quantum computer. First, we will find the quantum circuit for two spins and later extend it to $N$ spins with nearest neighbour interactions in one dimension. Each spin can be mapped to a qubit and the evolution of the spin system can be mapped to a quantum circuit. Using Eq. (5) and (6)

	$\displaystyle\prod_{\alpha=x,y,z}$	$\displaystyle e^{iJ_{\alpha}t(\sigma_{1}^{\alpha}\otimes\sigma_{2}^{\alpha})/\hbar}=$		(9)
		$\displaystyle\begin{pmatrix}e^{i\theta_{z}}\cos(\gamma)&0&0&ie^{i\theta_{z}}\sin(\gamma)\\ 0&e^{-i\theta_{z}}\cos(\delta)&ie^{-i\theta_{z}}\sin(\delta)&0\\ 0&ie^{-i\theta_{z}}\sin(\delta)&e^{-i\theta_{z}}\cos(\delta)&0\\ ie^{i\theta_{z}}\sin(\gamma)&0&0&e^{i\theta_{z}}\cos(\gamma)\end{pmatrix}$		(9)

where, $\gamma=\theta_{x}-\theta_{y}$ and $\delta=\theta_{x}+\theta_{y}$ . Optimal circuit for the above matrix is

\displaystyle\prod_{\alpha}

\displaystyle e^{iJ_{\alpha}t(\sigma_{1}^{\alpha}\otimes\sigma_{2}^{\alpha})/\hbar}=

(10)

\displaystyle\Qcircuit@C=0.5em@R=.7em@!R{&\ctrl{1}}

\displaystyle\gate{R_{x}(2\theta_{x})}

\displaystyle\gate{H}

\displaystyle\ctrl{1}

\displaystyle\gate{S}

\displaystyle\gate{H}

\displaystyle\ctrl{1}

\displaystyle\gate{R_{x}(-\pi/2)}

\displaystyle\qw

{\displaystyle\qw}

The evolution operator for any time step can be represented using the above circuit. In addition, it is also a constant depth circuit for each time step since the number of one- and two-qubit gates does not increase with time step. The quantum circuit for a spin chain with more than two spins in one dimension can be derived using Eq. (10) and the Trotter decomposition.

There exists two commuting families of operators as shown in Fig. 2 as orange and blue two qubit gates, respectively.

The accuracy of the simulation for a given time depends on the Trotter step ( $t/n$ ). As a consequence of Trotterization, the quantum circuit for time evolution grows linearly with time step. Fig. 2 shows the quantum circuit for a given time $t$ using $n$ Trotter steps. Each Trotter step is composed of a bilayer of two-qubit gates. The first layer acts on the first two qubits, followed by the third and then the fourth qubits and so on. Orange rectangles in Fig. 2 represent the first layer. The second layer of two qubit gates starts from second qubit and acts on the next two qubits. Blue rectangles in Fig. 2 represent the second layer. Both orange and blue rectangles combine to form an alternative layer, covering all possible nearest neighbour interactions.

IV Circuit Compression Using the Yang-Baxter Equation

In section III, we showed the generalized circuit for time evolution of $N$ spins in one dimension. In this section we will utilize the YBE to simplify the generalized quantum circuit for arbitrary time step. First, we will show the existence of a unique reflection symmetry for a quantum circuit composed of alternative layers. Next, we will show the merge identity for two-qubit gates. We will also show how reflection symmetry combined with the merge identity allows for the compression of a quantum circuit of any length to a finite depth.

IV.1 Reflection Symmetry and Merge Identity

The evolution operator for the Heisenberg Hamiltonian on two qubits is given by Eq. ( 9). When there are two evolution operators with different parameters on the same two qubits, they can be merged and represented via a single operator as shown below

\mathcal{R}^{ij}(\theta_{x}^{1},\theta_{y}^{1},\theta_{z}^{1}).\mathcal{R}^{ij}(\theta_{x}^{2},\theta_{y}^{2},\theta_{z}^{2})=\mathcal{R}^{ij}(\theta_{x}^{1}+\theta_{x}^{2},\theta_{y}^{1}+\theta_{y}^{2},\theta_{z}^{1}+\theta_{z}^{2})

(11)

In rest of this article we shall call Eq. (11) as the merge identity. Diagrammatic representation of the YBE is shown in the top panel of Fig. 3, where $\mathcal{R}^{ij}$ represents the operator acting on $i,j$ qubits. By using this symmetry repeatedly one can prove the existence of reflection symmetry in $n$ qubits composed of $n/2$ alternative layers [83]. The bottom panel of Fig. 3 shows how reflection symmetry is achieved for four qubits using the YBE three times.

Reflection symmetry combined with the merge identity allows for the compression of $N$ alternative layers of gates to $N/2$ alternative layers for $N$ qubits. Fig. 4 shows the use of reflection symmetry combined with the merge identity for four qubits. A third alternative layer can be merged into the previous two layers. Therefore, any number of alternative layers can be compressed into two alternative layers.

IV.2 Algebraic Condition for Reflection Symmetry

In the previous section, we have shown that reflection symmetry is a sufficient condition for performing circuit compression. An interesting follow-up question would be “How would one know if reflection symmetry can be applied to a quantum dynamics simulation of a given Hamiltonian?” To answer this question, one needs to show if there exist algebraic relations of phases and rotations before and after the reflection. However, given a general time propagator such as (9), an exhaustive search and rigorous proof appears challenging. Due to the lack of rigorous algebraic relations, YBE-like relations can only be conjectured as shown in previous work (see for example Ref. 41), which makes the compression of the circuit a heuristic process. In this section, we show that some algebraic relations can be obtained rigorously for at least a few special cases of (9). In particular, we propose the following theorem.

Theorem I Given the time evolution operator that takes the following form

\displaystyle\mathcal{R}(\gamma,\delta)

\displaystyle=\begin{pmatrix}e^{i\delta}\cos(\gamma)&0&0&ie^{i\delta}\sin(\gamma)\\ 0&e^{-i\delta}\cos{\gamma}&ie^{-i\delta}\sin{\gamma}&0\\ 0&ie^{-i\delta}\sin{\gamma}&e^{-i\delta}\cos{\gamma}&0\\ ie^{i\delta}\sin(\gamma)&0&0&e^{i\delta}\cos(\gamma)\end{pmatrix},

(12)

then the following YBE holds

	$\displaystyle(\mathcal{R}(\gamma_{1},\delta_{1})\otimes\mathbb{1})(\mathbb{1}\otimes\mathcal{R}(\gamma_{2},\delta_{2}))(\mathcal{R}(\gamma_{3},\delta_{3})\otimes\mathbb{1})$
	$\displaystyle~{}~{}~{}~{}=(\mathbb{1}\otimes\mathcal{R}(\gamma_{4},\delta_{4}))(\mathcal{R}(\gamma_{5},\delta_{5})\otimes\mathbb{1})(\mathbb{1}\otimes\mathcal{R}(\gamma_{6},\delta_{6}))$		(13)

if and only if the following 16 relations between the $\gamma$ ’s and $\delta$ ’s are satisfied

$\displaystyle s_{\gamma_{2}}c_{\gamma_{1}-\gamma_{3}}c_{\delta_{1}-\delta_{3}}s_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}s_{\delta_{4}+\delta_{6}}c_{\delta_{5}}$	(14)
$\displaystyle c_{\gamma_{2}}c_{\gamma_{1}-\gamma_{3}}c_{\delta_{1}+\delta_{3}}s_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}s_{\delta_{4}+\delta_{6}}c_{\delta_{5}},$	(15)
$\displaystyle-s_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}s_{\delta_{1}-\delta_{3}}c_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}s_{\gamma_{4}-\gamma_{6}}c_{\delta_{4}+\delta_{6}}s_{\delta_{5}},$	(16)
$\displaystyle c_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}s_{\delta_{1}+\delta_{3}}c_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}c_{\gamma_{4}-\gamma_{6}}c_{\delta_{4}+\delta_{6}}s_{\delta_{5}},$	(17)
$\displaystyle s_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}c_{\delta_{1}-\delta_{3}}c_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}c_{\delta_{4}+\delta_{6}}c_{\delta_{5}},$	(18)
$\displaystyle c_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}c_{\delta_{1}+\delta_{3}}c_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}c_{\delta_{4}+\delta_{6}}c_{\delta_{5}},$	(19)
$\displaystyle-s_{\gamma_{2}}c_{\gamma_{1}-\gamma_{3}}s_{\delta_{1}-\delta_{3}}s_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}s_{\gamma_{4}-\gamma_{6}}s_{\delta_{4}+\delta_{6}}s_{\delta_{5}},$	(20)
$\displaystyle c_{\gamma_{2}}c_{\gamma_{1}-\gamma_{3}}s_{\delta_{1}+\delta_{3}}s_{\delta_{2}}$	$\displaystyle=c_{\gamma_{5}}c_{\gamma_{4}-\gamma_{6}}s_{\delta_{4}+\delta_{6}}s_{\delta_{5}},$	(21)
$\displaystyle s_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}c_{\delta_{1}-\delta_{3}}c_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}c_{\delta_{4}-\delta_{6}}c_{\delta_{5}},$	(22)
$\displaystyle c_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}c_{\delta_{1}+\delta_{3}}c_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}c_{\delta_{4}-\delta_{6}}c_{\delta_{5}},$	(23)
$\displaystyle s_{\gamma_{2}}s_{\gamma_{1}-\gamma_{3}}s_{\delta_{1}-\delta_{3}}s_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}s_{\gamma_{4}-\gamma_{6}}s_{\delta_{4}-\delta_{6}}s_{\delta_{5}},$	(24)
$\displaystyle-c_{\gamma_{2}}s_{\gamma_{1}-\gamma_{3}}s_{\delta_{1}+\delta_{3}}s_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}c_{\gamma_{4}-\gamma_{6}}s_{\delta_{4}-\delta_{6}}s_{\delta_{5}},$	(25)
$\displaystyle-s_{\gamma_{2}}s_{\gamma_{1}-\gamma_{3}}c_{\delta_{1}-\delta_{3}}s_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}s_{\delta_{4}-\delta_{6}}c_{\delta_{5}},$	(26)
$\displaystyle-c_{\gamma_{2}}s_{\gamma_{1}-\gamma_{3}}c_{\delta_{1}+\delta_{3}}s_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}s_{\delta_{4}-\delta_{6}}c_{\delta_{5}},$	(27)
$\displaystyle-s_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}s_{\delta_{1}-\delta_{3}}c_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}s_{\gamma_{4}-\gamma_{6}}c_{\delta_{4}-\delta_{6}}s_{\delta_{5}},$	(28)
$\displaystyle c_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}s_{\delta_{1}+\delta_{3}}c_{\delta_{2}}$	$\displaystyle=s_{\gamma_{5}}c_{\gamma_{4}-\gamma_{6}}c_{\delta_{4}-\delta_{6}}s_{\delta_{5}},$	(29)

where $s_{p}$ and $c_{p}$ denote $\sin{(p/2)}$ and $\cos{(p/2)}$ , respectively.

The proof of Theorem I is straightforward if we expand both sides of Eq. (13) and perform a term-by-term comparison. Based on this theorem, we show explicitly in Table. IV.2, where we show YBE analysis for six special Heisenberg Hamiltonians where reflection is accomplished algebraically. Here, the Hamiltonian operator $\hat{H}$ takes at most two terms from the set $\{H_{X},H_{Y},H_{Z}\}$ where

$\displaystyle H_{X}$	$\displaystyle=-\sum_{j=1}^{n-1}J_{x}\sigma^{x}_{j}\sigma^{x}_{j+1},$	(30)
$\displaystyle H_{Y}$	$\displaystyle=-\sum_{j=1}^{n-1}J_{y}\sigma^{y}_{j}\sigma^{y}_{j+1},$	(31)
$\displaystyle H_{Z}$	$\displaystyle=-\sum_{j=1}^{n-1}J_{z}\sigma^{z}_{j}\sigma^{z}_{j+1}.$	(32)

$H_{X}+H_{Y}$	$\tau=\gamma+\delta,\phi=\gamma-\delta$		( $\tau,\phi$ )-relation is analogous to Eqs. (14-29)
Hamiltonian	Time Propagator	Time Propagator Circuit	Necessary Conditions for YBE
$H_{X}$	$e^{-itH_{X}/\hbar}=\mathcal{R}(\gamma,\delta=0)$	\Qcircuit@C=0.5em @R=.7em @!R & \ctrl1 \gateR_x(-2 γ) \ctrl1 \qw \targ \qw \targ \qw	$s_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}=c_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}$ , $c_{\gamma_{2}}c_{\gamma_{1}+\gamma_{3}}=c_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}$
$H_{X}$	$e^{-itH_{X}/\hbar}=\mathcal{R}(\gamma,\delta=0)$		$s_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}=s_{\gamma_{5}}s_{\gamma_{4}+\gamma_{6}}$ , $c_{\gamma_{2}}s_{\gamma_{1}+\gamma_{3}}=s_{\gamma_{5}}c_{\gamma_{4}+\gamma_{6}}$
$H_{Y}$	$e^{-itH_{Y}/\hbar}=U_{1}\mathcal{R}(\gamma,\delta=0)U_{1}^{\dagger}$	\Qcircuit@C=0.5em @R=.7em @!R & \gateR_z(π/2) \ctrl1 \gateR_x(-2 γ) \ctrl1 \gateR_z(-π/2) \qw \gateR_z(π/2) \targ \qw \targ \gateR_z(-π/2) \qw	Same as above
$H_{Z}$	$e^{-itH_{Z}/\hbar}=\mathcal{R}(\gamma=0,\delta)$	\Qcircuit@C=0.5em @R=.7em @!R & \ctrl1 \qw \ctrl1 \qw \targ \gateR_z(-2 δ) \targ \qw	$c_{\delta_{1}+\delta_{3}}s_{\delta_{2}}=s_{\delta_{4}+\delta_{6}}c_{\delta_{5}}$ , $s_{\delta_{1}+\delta_{3}}c_{\delta_{2}}=c_{\delta_{4}+\delta_{6}}s_{\delta_{5}}$
$H_{Z}$	$e^{-itH_{Z}/\hbar}=\mathcal{R}(\gamma=0,\delta)$		$c_{\delta_{1}+\delta_{3}}c_{\delta_{2}}=c_{\delta_{4}+\delta_{6}}c_{\delta_{5}}$ , $s_{\delta_{1}+\delta_{3}}s_{\delta_{2}}=s_{\delta_{4}+\delta_{6}}s_{\delta_{5}}$
$H_{X}+H_{Y}$	$e^{-it(H_{X}+H_{Y})/\hbar}=U_{2}\mathcal{R}(\tau,\phi)U_{2}^{\dagger}$	\Qcircuit@C=0.5em @R=.7em @!R & \gateR_x(π/2) \ctrl1 \gateR_x(- 2τ) \ctrl1 \gateR_x(-π/2) \qw \gateR_x(π/2) \targ \gateR_z(- 2ϕ) \targ \gateR_x(-π/2) \qw	( $\tau,\phi$ )-relation is analogous to Eqs. (14-29)
$H_{X}+H_{Z}$	$e^{-it(H_{X}+H_{Z})/\hbar}=\mathcal{R}(\gamma,\delta)$	\Qcircuit@C=0.5em @R=.7em @!R & \ctrl1 \gateR_x(-2 γ) \ctrl1 \qw \targ \gateR_z(-2 δ) \targ \qw	Eqs. (14-29)
$H_{Y}+H_{Z}$	$e^{-it(H_{Y}+H_{Z})/\hbar}=U_{1}\mathcal{R}(\gamma,\delta)U_{1}^{\dagger}$	\Qcircuit@C=0.5em @R=.7em @!R & \gateR_z(π/2) \ctrl1 \gateR_x(-2 γ) \ctrl1 \gateR_z(-π/2) \qw \gateR_z(π/2) \targ \gateR_z(-2 δ \targ \gateR_z(-π/2) \qw	Eqs. (14-29)

$\displaystyle\tan([\gamma_{4}+\gamma_{6}]/2)$	$\displaystyle=\tan(\gamma_{2}/2)\frac{\cos([\delta_{1}-\delta_{3}]/2)}{\cos([\delta_{1}+\delta_{3}]/2)},$	(33)
$\displaystyle\tan([\gamma_{4}-\gamma_{6}]/2)$	$\displaystyle=-\tan(\gamma_{2}/2)\frac{\sin([\delta_{1}-\delta_{3}]/2)}{\sin([\delta_{1}+\delta_{3}]/2)},$	(34)
$\displaystyle\tan([\delta_{4}+\delta_{6}]/2)$	$\displaystyle=\tan(\delta_{2}/2)\frac{\cos([\gamma{1}-\gamma{3}]/2)}{\cos([\gamma{1}+\gamma{3}]/2)},$	(35)
$\displaystyle\tan([\delta_{4}-\delta_{6}]/2)$	$\displaystyle=\tan(\delta_{2}/2)\frac{\sin([\gamma{1}-\gamma{3}]/2)}{\sin([\gamma{1}+\gamma{3}]/2)},$	(36)
$\displaystyle\tan(\gamma_{5}/2)$	$\displaystyle=\tan([\gamma_{1}+\gamma_{3}]/2)\frac{\cos([\delta_{4}+\delta_{6}]/2)}{\cos([\delta_{4}-\delta_{6}]/2)},$	(37)
$\displaystyle\tan(\delta_{5}/2)$	$\displaystyle=-\tan([\delta_{1}+\delta_{3}]/2)\frac{\cos([\gamma_{4}+\gamma_{6}]/2)}{\cos([\gamma_{4}-\gamma_{6}]/2)}.$	(38)