Algebraic Bethe Circuits

The Algebraic Bethe Ansatz (ABA) is a highly successful analytical method used to exactly solve several physical models in both statistical mechanics and condensed-matter physics. Here we bring the ABA into unitary form, for its direct implementation on a quantum computer. This is achieved by dis-tilling the non-unitary R matrices that make up the ABA into unitaries using the QR decomposition. Our algorithm is deterministic and works for both real and complex roots of the Bethe equations. We illustrate our method on the spin-12 XX and XXZ models. We show that using this approach one can eﬃciently prepare eigenstates of the XX model on a quantum computer with quantum resources that match previous state-of-the-art approaches. We run small-scale error-mitigated implementations on the IBM quantum computers, including the preparation of the ground state for the XX and XXZ models on 4 sites. Finally, we derive a new form of the Yang-Baxter equation using unitary matrices, and also verify it on a quantum computer.


Introduction
One of the most widely-recognized applications of quantum computing is the efficient simulation of many-body quantum systems [1]. Simulating such systems using classical devices generally requires computational resources that scale exponentially with the size of the system. Quantum computers on the other hand are naturally suited to this task, being themselves quantum systems and hence overcoming the classical exponential scaling. In this context, preparing Bethe Ansatz (BA) eigenstates on a quantum computer is attracting increasing attention [2][3][4][5].
The BA is an extremely successful classical method for exactly solving one-dimensional (1D) quantum models, e.g. the Heisenberg, Hubbard or Kondo models [6][7][8][9]. It reduces the difficult problem of diagonalizing the Hamiltonian to finding the solutions of a * The first two authors contributed equally. Figure 1: A Bethe ansatz eigenstate |ΨM for M magnons on N sites, converted into a quantum circuit of N qubits by computing iterative QR decompositions and removing the ancilla qubits, as explained in the main text. The Rj are the tensors defined in (6) that satisfy the YB equation, and the P k are unitary matrices of dimension 2 n+1 × 2 n+1 with n = min(k, M ). These matrices need to be compiled to quantum gates before implementation on a quantum computer.
set of algebraic equations. In many cases, it is possible to numerically solve these equations, which in turn allows one to calculate the eigenvalues and the eigenvectors of the system of interest. These quantities are computed differently depending on whether the coordinate Bethe Ansatz (CBA) or the Algebraic Bethe Ansatz (ABA) is used [9]. In both cases, the eigenstates are represented as complex mathematical expressions. As a result, this method does not allow direct access to some physical quantities, such as high-order and long-range correlation functions, which have proved challenging to compute both analytically and numerically [10, and references therein]. This motivates the construction of such states directly on a quantum computer. Once the eigenstates are experimentally available, all the correlation functions can be readily computed from measurements. Furthermore, these states can be used to initialise other quantum algorithms or benchmark quantum hardware.
To this end, a quantum algorithm for the preparation of the eigenstates of the spin-1 2 XXZ model in 1D was introduced in Refs. [3,4], based on the CBA. This algorithm is probabilistic and works for real-valued solutions of the Bethe equations. Its circuit depth is polynomial in both the system size and the number of magnons, or down spins. However, it was recently shown [5] that the success probability of obtaining the desired eigenstate upon measurement of the ancillary qubits decreases super-exponentially with the number of magnons for large chains. The use of a variational approach for the preparation of integrablesystem eigenstates, which does not rely on knowledge coming from the BA, has also been considered [11][12][13][14]. Variational quantum algorithms are known to suffer from exponentially-vanishing gradients [15,16]. This problem is certain to appear when the number of magnons scales with the system size [17]. Moreover, even when the number of magnons is fixed this issue is likely to persist for large chains due to the effect of noise [18].
Another approach could be to prepare eigenstates of interest using Matrix Product State (MPS) optimization [19] and then use further optimization techniques to map the classically tractable state to a quantum circuit. However, calculating an accurate classical representation depends on the structure of entanglement of the systems and will not be trivial in all phases or the entire spectrum of states [20]. Furthermore, once a classical representation of the state of interest is calculated constructing a faithful quantum circuit to represent it in a quantum computer may also be difficult and is an active field of research [21,22].
In this paper, we present a quantum algorithm based on the ABA for the preparation of BA eigenstates. Contrary to [3], it also works for complex solutions of the Bethe equations. The main difficulty encountered when directly trying to convert the ABA into a quantum circuit is that the matrices R are not unitary. Besides, the translation of the ABA into a circuit requires the use of ancillas which need to be projected onto the |0 . . . 0 state at the end of the computation. This would yield a probabilistic algorithm should a direct translation be employed. We circumvent both of these issues by iteratively com-puting the QR decomposition [23] of the non-unitary matrices appearing in the ABA. This then allows us to obtain the desired eigenstate as the output of a quantum circuit with no ancillary qubits, which we refer to as an Algebraic Bethe Circuit (ABC). The process is sketched in Fig. 1. We note that quantum circuits based on unitary R matrices were explored to compute infinite temperature correlation functions in [24].
We focus on the paradigmatic 1D antiferromagnetic spin-1 2 XXZ model with periodic boundary conditions for the concrete analysis. The complexity of calculating the ABC unitaries scales linearly in the number of sites (qubits), but exponentially in general with the number of magnons. There is the additional complexity of compiling the calculated unitaries, which may also scale exponentially. This means that our algorithm is in general only applicable to a small number of magnons. In spite of that, it could allow for the preparation of states on near-term quantum hardware that have been challenging before. Besides, we find that its application on the XX model is efficient. This model describes free spinless fermions via a Jordan-Wigner transformation. Here, the ABCs match the performance of state-of-the-art algorithms for the preparation of fermionic states [25][26][27], both in the number of gates necessary and the circuit depth.
From a theoretical standpoint, the ABCs offer an alternative approach towards finding exact circuits for integrable quantum many-body systems [28,29]. Along these lines, we derive a novel version of the Yang-Baxter equation in terms of unitary matrices which can be tested on quantum hardware. Interestingly, the BA can be interpreted as a MPS [30][31][32]. Our algorithm for the distillation of unitaries from the ABA appears closely related to the transformation of an MPS into canonical form [33], which in turn has a direct translation to a quantum circuit. Hence our method to obtain ABCs should also prove relevant for the circuit implementation of general MPS with low bond dimension. Let us mention that several works have considered the implementation of tensor-network states on a quantum computer [22,34,35,[35][36][37][38].
All current quantum computers suffer from significant hardware noise. Therefore, error mitigation will play an essential role in all near-term quantum algorithms and simulations. As such, we explore the performance of several state-of-the-art error mitigation methods on mitigating observables produced by our algorithm.
The paper is structured as follows. In Section 2 we briefly introduce the ABA. Section 3 contains a detailed explanation of its transformation into a quantum circuit, together with a preliminary discussion on the decomposition of the ABC unitaries in terms of elementary quantum gates. The unitary version of the YB equation is discussed in Section 4. Section 5 contains small-scale error-mitigated implementations of the ABCs on quantum hardware. In particular, we prepare plane wave states on 8 sites and implement the ground state of the XX and XXZ models on 4 sites using the IBM cloud quantum computers. Section 5 also includes a test of the YB equation on quantum hardware. Finally, we conclude in Section 6 with a discussion of our results. Several technical details are confined to Appendices.

The Algebraic Bethe Ansatz
The Algebraic Bethe Ansatz [7-9] is a powerful classical technique to solve one-dimensional quantum integrable vertex models. These models fulfill the socalled Yang-Baxter (YB) equation [39], and are characterised by an extensive set of conserved quantities. They have been greatly studied over the last century in the context of quantum many-body physics [40,41].
The ABA can be used to calculate the exact eigenspectrum of a large class of Hamiltonians. Among these models, a prominent example is the 1D spin-1 2 anti-ferromagnetic XXZ Hamiltonian on N sites with periodic boundary conditions, where {σ X j , σ Y j , σ Z j } are the Pauli matrices acting on the j-th spin, and σ N +1 ≡ σ 1 . The parameter ∆ introduces an anisotropy in the chain. When ∆ = 1 we recover the isotropic anti-ferromagnetic Heisenberg spin chain or XXX model, which has SU (2) symmetry. For other values of ∆, only a U (1) symmetry is present. A set of algebraic equations needs to be solved in order to construct the eigenstates of the Hamiltonian. These are the celebrated Bethe equations, whose solutions or roots are referred to as rapidities which we denote as {λ j } M j=1 . The number of Bethe roots M describes the number of magnons or spin-down waves composing the state.
We shall consider the XXZ model with anisotropy in the range ∆ ∈ (−1, 1]. The corresponding Bethe equations are given by where cos γ = ∆, with γ ∈ [0, π). The ground state of a chain with an even number of sites is built out of M = N/2 magnons whose rapidities are real. In general, excited eigenstates however contain complex rapidities, that must come in conjugate pairs to guarantee that the energy is real-valued. The Bethe equations are typically solved by numerical methods [42].
The BA eigenstates can be represented as a product of operators acting on the vacuum state |vac as In the XXZ model, the |vac state is the product state with all spins up, i.e. |0 . . . 0 , and B(λ) is an operator that creates one magnon. This operator can be represented as the contraction of a network of fourindex tensors (see Appendix B for a more detailed explanation) .
with the last equality making clear the matrix interpretation of the tensors appearing in (4). The input and output indices ai and jb correspond respectively to the columns and rows of the R matrix. The R matrix of the XXZ model is with ρ a complex number and the parameters s 1 , s 2 satisfying 1 + s 2 2 − s 2 1 = 2s 2 ∆ .
In terms of the rapidity, this equation is solved by The second parameter has a direct physical interpretation as the magnon quasi-momentum, s 2 = e ip . It is important to note that the R matrix is excitation preserving as a consequence of the U (1) symmetry of the Hamiltonian (1).
3 From ABA to ABC

Detailed method
Our aim is to construct quantum circuits based on the ABA, for the direct preparation of eigenstates of integrable vertex models on quantum hardware. We will call these circuits Algebraic Bethe Circuits (ABCs).
An ABA eigenstate for M magnons on N sites is shown in Fig. 1, where it is transformed into a suggestive form with the apparent structure of a quantum circuit. The basic cell R T that repeats itself throughout this circuit is , (9) where R j ≡ R(λ j ). The problem encountered when directly trying to transform this cell into a quantum gate is that the matrices R are in general not unitary (see Appendix A).
To stress the crucial difference between unitary and non-unitary matrices, we shall use rounded-corner rectangles for the latter. The complete ABA network can be recast in terms of the R T cells as The M rightmost qubits, both at the input and output, are in the fixed state |0 . They can be considered ancillary qubits. Keeping them in the final quantum circuit we are seeking would result in a probabilistic algorithm. A severe reduction of the success probability with increasing M is to be expected due to the curse of dimensionality. The probabilistic quantum algorithm based on the CBA proposed in [3] was recently shown to suffer from a similar problem [5].
Both the conversion of (10) into a quantum circuit and the removal of the ancillary qubits can be addressed by utilising the QR decomposition as our main tool. To be more precise, any m × n matrix can be written as the product of two matrices Q·R. When m ≤ n, Q is a m × m unitary and R an m × n rectangular matrix with vanishing entries below the main diagonal. When m > n, Q is a m × n isometry (i.e. Q † Q = 1 n ) and R an n × n upper triangular matrix. We note that the QR decomposition was previously used in the derivation of efficient quantum circuits in [25][26][27]43].
We start our protocol at the top-rightmost basic cell. All the information in this cell can be encoded in a 2 × 2 M matrix defined by . (11) Working directly with G 0 allows us to eliminate the rightmost ancillary qubit. The matrix G 0 can then be absorbed into the second R T cell. This defines a 4 × 2 M matrix which renders the second rightmost qubit unnecessary. We apply now the QR decomposition to this matrix, obtaining We have distilled the first gate of the deterministic quantum circuit for the construction of Bethe eigenstates, the two-qubit unitary P 1 . The non-unitary remainder G 1 is then absorbed into the next basic cell and the QR decomposition is computed again for the new non-unitary matrix. As before, one more ancillary qubit is eliminated. This process is iterated. At each step a unitary gate P k acting on k + 1 qubits is obtained. After M − 1 iterations all ancillary qubits have been removed and the M top-right basic cells have been substituted by the circuit . (13) For k ≥ M each new step is described by the recursion relation or equivalently The LHS defines a 2 M+1 × 2 M matrix. The matrix Q resulting from its QR decomposition is in this case an isometry, which determines P k |0 . This information can be completed at our best convenience to define the M + 1 qubit gate P k that is to be implemented on the quantum circuit. In order to solve (15), we multiply both sides with the Hermitian conjugate to obtain a recursion relation involving only the upper triangular matrices G k , i.e.
Once the solution G k has been calculated, we extract the P k from For k < M a similar relation holds Recall that G k<M−1 are not square matrices and hence do not have an inverse. However, there exists a non-square upper triangular matrix H k satisfying G k H k = 1, which is thus sufficient for our purposes.
Although the dimensions of G k vary in the first iterations, the product G † k G k always defines a 2 M × 2 M matrix. This renders the previous equation applicable to all iterations.
The protocol ends at the bottom-left cell The correct normalization of the output Bethe state |Ψ M from the ABA requires one to suitably adjust the global factor ρ of the R matrix (6). If this is done, g will be a trivial phase. Alternatively, we might ignore the difficult issue of finding the appropriate ρ for the desired eigenstate and just set ρ = 1. The output state from (10) will not be properly normalized, but the problem is simply solved by discarding the global factor g. We adopt this convention in the following. The previous relation (19) only determines the action of P N −1 on the state |1 . . . 10 . This last unitary can be otherwise chosen freely. We have thus obtained a deterministic quantum circuit for the preparation of BA eigenstates on N sites requiring only N qubits. The upper triangular structure of G k is key to the feasibility of our protocol. It implies the reduction of (16) to a nested system of equations. Substituting the solution of one equation into the next, we only need to care for one entry of G k at a time. Moreover the inverse of a triangular matrix, required to obtain the unitaries P k , is not expensive to calculate. The bottleneck of the procedure is the exponential growth in the number of equations with the number of magnons. Nevertheless, we have found that this is not prohibitive to explore cases of interest.
Two important comments on the QR decomposition should be added. The matrices Q and R are not fully determined but enjoy the gauge freedom where D is an arbitrary diagonal matrix containing complex phases. This freedom will be crucial below, when compiling the unitaries to elementary quantum gates. Additionally, the QR decomposition is compatible with the U (1) symmetry of the XXZ model. Hence, all matrices G k and P k will conserve individually the number of excitations. We end with a final remark. It has been shown that the BA can be interpreted in the language of tensor networks as a Matrix Product State (MPS) [30][31][32]. Our algorithm for the distillation of unitaries from the R T matrices is closely related to the transformation of an MPS into canonical form [33]. The tensors of a canonical MPS define isometries, and thus have a direct translation to a quantum circuit. The main difference with our approach lies in the elimination of the ancillary qubits, which renders the circuit deterministic. Hence the ABCs could prove relevant, not only for the preparation of eigenstates of integrable spin chains, but also for general MPS of low bond dimension. We note that the implementation of tensor-network states on a quantum computer has been addressed in [22,34,35,[35][36][37][38].

One-magnon solution
We illustrate our method by explicitly deriving the circuit that constructs one-magnon states for the XXZ model. In order to study the properties of the P k gates, we find it convenient to use the magnon quasimomentum p instead of the rapidity λ as the basic variable. Indeed, the R matrix (6) has a simple expression in terms of s 1 and s 2 = e ip . Besides, the integrability constraint (7) turns the parameter s 1 into a function of the quasi-momentum and the anisotropy.
The basic cell constructing one-magnon solutions is composed of a single R matrix. In this case all the G k matrices are 2 × 2 upper triangular. The preservation of U (1) symmetry further forces them to be diagonal. We choose the following ansatz The initial matrix G 0 defined in (11) is of this form with c 0 = 1. For k > 0, the gauge freedom of the QR decomposition allows us to set the c k parameters to be real and positive. Equation (16) translates then into the recursion relation for k ≥ 1. Using the initial condition c 0 = 1, we arrive at the simple solution where we have used that the quasi-momentum of the one-magnon solutions must be real, implying that s 2 is a phase. Substituting into (17) and completing the matrix, we obtain the two-qubit quantum gates written with the convention We have ordered the input and output qubits from right to left instead of the customary left to right order, since this is better suited to the structure of the ABC. Indeed, with this convention the first two columns of (24) describe P k |0 . One-magnon solutions are just plane waves. The anisotropy parametrizes the strength of interactions among magnons and it hence should have no influence on these configurations. Consequently, the unitaries P k are independent of ∆. This is in contrast with the ABA circuit network (10), where every element depends on both s 1 and s 2 , or equivalently, on ∆ and p. It is only the output of the complete circuit that builds the one-magnon states which is independent of the anisotropy. In our scheme, the local dependence on ∆ is confined to the unphysical matrices G k . Our construction offers an economic way of implementing plane waves on quantum hardware, with just nearestneighbor connectivity among qubits. The special case p = 0 reproduces the optimal algorithm for obtaining the W state [44].

The XX model
The XX model, obtained when ∆ = 0, is the simplest member of the XXZ family. It describes a spinless free fermion system via the Jordan-Wigner transformation. This was used in [28] to construct an efficient quantum circuit to prepare its eigenstates with a number of gates that is quadratic in the number of sites N . We will show that the ABC requires O(N M ) gates.
We search for a decomposition of the unitaries P k in terms of two-qubit quantum gates. In particular, we will use the phased fSim (F ) gate [25] as the basic building block of our circuit, defined as This matrix is a U (1)-symmetry-preserving generalization of the one-magnon unitaries (24). It can be implemented with five R X gates, five R Z gates and two CNOTs. We note that this unitary is a so called matchgate [45]. Circuits consisting of matchgates are classically simulable and correspond to systems of free fermions. We use an ansatz in which the F gates reproduce the same contraction structure as the R matrices in the basic R T cell, i.e.
. (27) In Appendix D we present the complete solution for two and three magnons in the XX model, including the closed-form of the matrices P k . Furthermore, we verify the previous ansatz and analytically derive the parameters of the F gates. We have checked numerically that (27) holds up to six magnons. It also holds for the unitaries P k with k < M , which act on a reduced number of qubits. Based on this evidence, we conjecture that it is valid in general for the XX model. In spite of the similar decompositions of the matrices P k and R T , there is an important difference between them. Each R j matrix in R T depends on a single magnon quasi-momentum p j . On the contrary, the two-qubit gates F j,k contain information from all magnons down to position j, as seen in (73)-(75). Explicitly, the gate F M,k is a function of a single quasi-momentum p M , while F 1,k involves all of them.
The phase freedom shown in (20) is essential for (27) to hold. It is simple to see that the onespin-down sector of P k fixes completely this ansatz. The free phases affecting other symmetry sectors have then to be chosen appropriately for the so derived F j,k to describe the complete P k . This choice has been implemented in the two and three magnon unitaries presented in Appendix D. Since G k and P k preserve the U (1) symmetry, the one-spin-down sector can be determined independently from all others from the recursion relation (16). This implies that only M (M +1) 2 classical equations need to be solved. Hence our algorithm for the XX model is efficient both in its classical and quantum parts.
The total number of two-qubit gates for an ABC creating M magnons on N sites is N M − M (M +1) 2 .
All gates act upon nearest-neighboring qubits, which means that it is possible to implement it on quantum hardware with simple 1D connectivity. The total depth of the circuit is N + M − 1, which represents a factor O(log(N )) improvement over the method presented in Ref. [28]. Furthermore, it would be interesting to explore whether the depth of the ABC can be further reduced by considering additional connectivity [44].
It is worth noting that the circuit that we have obtained for the XX model is closely related to those considered in Refs. [25][26][27] for free fermion models, that are based on preparing Slater determinants. These circuits have the same behavior in the depth and number of gates as the ABCs. This is to be expected since this model is non interacting and the eigenstates are themselves Slater determinants. There are sure to be interesting ties between the two approaches that may prove a fruitful research direction.

The XXZ model
The XXZ chain with non-vanishing anisotropy is an interacting model and thus computationally more challenging. Therefore, we expect the unitaries P k making up our ABC in the XXZ model to require a number of gates that grows exponentially with the number of magnons. Moreover, the number of equations that need to be solved in order to find the full matrices also scales exponentially with M . However, for sufficiently small numbers of magnons we find that this is not prohibitive. Indeed, we can numerically obtain the matrices P k up to M = 12 with modest computational resources. It should be noted that whether it is strictly necessary or not to calculate the entire matrices in order to find a quantum circuit that implements them is a matter of investigation. Likewise, whether an efficient quantum circuit for increasing M may exist cannot be completely discarded based on our results.
For the decomposition of the unitaries P k we define theF gates identically to the F gates (26) but include an additional phase e iγ in the last diagonal element of the matrix. We note that the addition of this phase leads to the gates no longer being matchgates. Therefore, including this phase means the circuit can no longer be mapped to free fermions. This gate can be implemented with ten one-qubit gates and three CNOTs [46]. When M = 2, we find that it is possible to decompose P k with . (28) The error in the compilation is measured using achieving an order 10 −10 precision. This renders the preparation of two-magnon configurations for small systems accessible to current quantum hardware. For a larger number of magnons there is more freedom in how to distribute theF gates. As a first exploratory ansatz, we use the repetition of the structure in (27) in a layer-wise fashion, allowing for γ phases in all gates.
We have verified up to 5 magnons that the compilation error associated with this ansatz does not depend on the matrix P k for k ≥ M . Therefore, we have focused on P M , i.e. on the unitary obtained right after the elimination of the ABA ancillary qubits. Its numerical decomposition is addressed for several values of the anisotropy and M = 3, 4, 5, as shown in Fig. 2. Our results indicate that a very small compilation error is obtained with around 8, 18, 55 layers respectively, making up a total of 24, 72, 275F gates. The number of layers therefore appears to grow exponentially with the number of magnons for this ansatz. The compilation was achieved using the Berkeley Quantum Synthesis Toolkit (BQSKit) [47], with a slight modification introduced to account for the fact that we are actually compiling P M |0 instead of the full unitary matrix. It should be stressed that these numerical results constitute a very preliminary study of the decomposition of the unitaries in the ABCs.
As a complementary direction, it is important to analyze the general properties of the unitaries P k . With this aim we have derived the complete analytical solution for M = 2 in the XXZ model (see Appendix C). We have seen that the first columns of P k reproduce the one-magnon solution (24). Also, we have checked that the same property holds for three magnons, and conjecture that it holds in general. With this assumption, the solutions with M < M magnons turn out to be contained in those for M magnons. This property is evident for R T and the P k 's from the XX model due to their single-layer structure of two-qubit matrices preserving the number of excitations. Since this simple structure is lost in the XXZ unitaries, we find remarkable that the hierarchy of solutions is preserved.
The simple structure of R T and the XX model P k 's has another important consequence. It forces many entries of these matrices to vanish in spite of being compatible with the U (1) symmetry. Starting from M = 3, some of these entries become non-zero when the anisotropy does not vanish. This is at the core of the difficulty encountered when decomposing the unitaries of the XXZ model. A better understanding on the emergence of these entries should be crucial in searching for an optimal decomposition. Furthermore, the dependence on the index k is an essential aspect of the ABCs, which deserves study both because of its theoretical and practical implications. A detailed investigation along these lines is left for future work.
Finally, it is interesting to note that the distinction between the XX model having a polynomial-depth circuit in both the number of qubits and the number of magnons, and the XXZ model likely having a circuit depth that scales exponentially with the number of magnons of the eigenstate, is matched by the dimension of the dynamical Lie algebras of the corresponding Hamiltonian generators. For the XX model, the dimension of this algebra is polynomial in the number of sites irrespective of the excitation subspace [48]. In contrast, for the XXZ model the dimension of the subspace algebra is exponential when the number of magnons scales linearly with the number of sites [17]. Whether there exists a direct causal relation between these observations is left for now as an open question to be explored in future work.

Unitary form of the Yang-Baxter equation
The R matrices of an integrable quantum system satisfy the celebrated YB equation where λ and µ are the rapidity parameters introduced in Section 2. This equation is the statement that any process involving N sites can be factorized into twobody pieces. It guarantees that the BA state given in (3) is independent of the order of the parameters λ 1 , . . . , λ M . Let σ be an arbitrary permutation of the M rapidities. There exists then a 2 M × 2 M tensor R σ satisfying , (31) where R σ T is defined by (9) with the permuted rapidities. When substituting it into the network (10), the independence on the ordering of the rapidities becomes manifest. The simplest case of this relation is provided by the YB equation itself, whose graphic representation is . (32) Here σ permutes λ and µ, and R σ = R(λ − µ). For solutions containing more than two magnons, R σ consists of a product of R matrices depending on differences of rapidities. Although in general this product is not unique, the basic equation (32) ensures that all choices lead to the same matrix R σ (see Appendix F).
We now analyze the implications of the YB equation on the quantum gates defining the circuit version of the ABA. Relation (31) translates straightforwardly into . (33) The exchange matrices M k are obtained by dressing R σ with the non-unitary pieces left by the QR decomposition that brings the ABA into unitary form with G σ k depending on the permuted rapidities. Moreover, (31) also implies that the product G σ k R σ satisfies the same recursion relations (16) that define G k . Both sets of solutions are related by which is proven by induction, using that G σ 0 R σ = G 0 . It immediately follows that the exchange matrices M k,σ are unitary. Hence (33) is the unitary version of the YB equation. It describes how the basic cell of the Bethe circuit changes under a permutation of the rapidity parameters. This reformulation of the YB equation allows for its test on quantum hardware with no restriction on the rapidities (see below). It represents another application of the YB equation, and adds to those found in exactly-solvable models in Statistical Mechanics [39] and the factorized S-matrices in relativistic quantum field theory models [49].
The exchange matrices for M = 2 are calculated from the direct relation with the R matrix . (36) We have dropped the sub-index σ in M k since in this case there is only one possible permutation, which interchanges λ and µ. The analytical expression is presented in Appendix G. Although all elements on the right-hand side of (36) depend on the anisotropy, their product does not. This implies that M k is the same function of the magnon quasi-momenta for all members of the XXZ family. This shows that the matrices (36) cannot be interpreted as unitary versions of the R matrix. We stress that the exchange matrices depend on the anisotropy for M ≥ 3. These matrices are interesting objects in their own right, deserving further study.

Numerical and Experimental Results
We performed numerical simulations to verify our theoretical results, using the open-source library Qibo [50,51]. We numerically solved the Bethe equations (2) and simulated the unitary circuits at the level of the P k matrices for several excited states of the XXZ model up to 24 sites and 12 magnons. We compared the resulting eigenstates with the simulated ABA states, finding a perfect agreement between the two. The circuits were simulated in double precision using the qibojit backend [52] on multi-threading CPU. They were directly obtained by computing the QR decompositions that are required to convert the ABA into a deterministic quantum circuit. The programs to reproduce such simulations can be found at [53].

Plane waves on quantum hardware
We also implemented the circuit construction for plane wave states for a system size N = 8 on the quantum computer IBM_Montreal. We explored the twopoint correlators from the prepared state (see Fig. 3).
As previously shown, the 2 qubit P k unitary gates are simply phased fSim gates when M = 1. These gates were decomposed into the IBM native gate set to be implemented on the hardware. The circuits used to prepare these states have total depth 65. Each unitary gate involved 2 CNOT gates, leading to 2(N − 1) CNOT gates in total. Current devices suffer from significant hardware noise. In order to obtain the best possible results it is necessary to use error mitigation which focuses on reducing the impact of noise. There are many current error mitigation approaches and unifications thereof, each with their respective advantages and disadvantages [54,55]. Exploring how error mitigation performs for a task of interest on real quantum computers give insight onto their performance in realistic scenarios.
In this work we implemented three techniques: zero-noise extrapolation (ZNE) [56], Clifford data regression (CDR) [57] and variable noise Clifford data regression (vnCDR) [58]. We used the open source software package Mitiq [59] to execute these methods. For more details regarding the implementation of these techniques we refer the reader to Appendix H.
We benchmark the performance of each method by calculating a weighted error, defined as: where O j is some observable of interest and O j exp is the estimated value for that observable obtained experimentally with or without error mitigation. The mean is taken over the different qubit positions. This definition prevents the error from diverging if O j exact approaches zero while allowing the errors corresponding to several observables to be compared.
To simplify the comparison between mitigation techniques we average the above error metric across all the different observables. We expect this metric to reflect the overall performance of the mitigation methods while also taking into account the magnitude of the observables to mitigate. Furthermore, to simplify the presentation we omit the CDR mitigated results from our plots, we see that in general the performance is significantly better than ZNE but worse than vnCDR. Error mitigation significantly improves the results obtained from the real device (see Fig. 3). As the observables become less local and the circuit depth increases, the quality of the raw and error mitigated results tends to decrease. This is particularly noticeable in the results obtained for the σ Z 1 σ Z j correlators. Clearly both vnCDR and ZNE tend to perform best in shallower, less noisy circuits although they still improve results for the deepest circuits explored here.

Two-magnon states on quantum hardware
In addition, we implemented the circuits for two magnon states on current quantum hardware. We constructed the ground state of the XXZ model for 4 sites and an excited state of the XX model for 5 sites. We constructed the ground state of the XX model for 4 sites in Appendix E.
The circuit for the ground state of the XXZ model consists of 5 F gates and 2F gates, structured as  follows: . (38) We only need to determine the action of the last unitary P N −1 on the state |011 and this can be achieved with a single layer of F gates.
Once the F andF gates have been decomposed into the IBM native gate set the circuit to prepare these states have depths 57 and 49 and involve 16 and 14 CNOT gates respectively.
We evaluated the two-point correlators for these states across the chain, and mitigated the effect of hardware noise with ZNE and vnCDR (see Fig. 4). The mitigation methods improve upon the raw data. For two magnon states we explored smaller system sizes due to the increased scaling of depth with the number of qubits. We find that good agreement can be obtained between the exact and error mitigated observables, with vnCDR reducing the effect of noise the most.
Overall, these experiments show a proof of principle implementation of our approach for low numbers of magnons. Furthermore, they highlight the utility of error mitigation. In particular they show further evidence that learning based error mitigation is practically useful in reducing the effects of hardware noise. For larger scale implementations a combination of noise reduction and error mitigation techniques will be needed, which presents an exciting challenge for fu-  Table 1: Verification of the unitarised YB equation using the IBM_Cairo cloud quantum computer. The left column shows the initial state fed into (33), and the right column shows the fidelity between the output states at both sides of this equation.

Yang-Baxter equation on quantum hardware
The significance of the YB equation has motivated its experimental implementation. A two-dimensional version of the YB equation has been verified optically [60] and by nuclear magnetic resonance [61]. We have also verified the unitarised YB equation (33) on the cloud quantum computer IBM_Cairo, for the case M = 2 magnons in the XX model, given by We used state tomography to compute the density matrices associated with the right (ρ r ) and left (ρ l ) output states at both sides of (39) (see Table 1). We computed these density matrices for each of the four possible initial input states, and we did so for the P 2 gate with Bethe roots λ 1 = −1/ √ 3 and λ 2 = 1/ √ 3. In order to compare ρ r and ρ l we determined the fidelity, given by This heralds the first implementation of the YB equation on a superconducting quantum computer.

Discussion
In this work, we have introduced a method to exactly prepare eigenstates of quantum integrable vertex models on programmable digital quantum computers. Our approach relies on using the QR decomposition as the main tool to bring the ABA to unitary form. In contrast to previous proposals, our method works for both real and complex roots of the Bethe equations and is deterministic. Both the circuit depth and gate complexity of our approach scale linearly with the number of qubits. However, we expect an exponential scaling in the number of magnons in general. This could affect the classical preprocessing, the required compilation step, and the circuit depth. Despite this, we find that with modest classical computational resources one can obtain a unitary circuit representation for interesting states. For the XX model, which can be mapped to free fermions, we find an efficient gate decomposition with polynomial classical effort. In particular, our approach produces quantum circuits that match the state-of-theart O(N ) depth of [25][26][27]62].
Our algorithm opens up the possibility to prepare highly non-trivial ABA eigenstates on quantum computers. Foreseeable applications include using these states to study Hamiltonian quenches that may be inaccessible to classical methods. More generally, these states could be used as inputs to other quantum algorithms. For instance, it would be interesting to explore if they can be used to initialize variational quantum algorithms. Using such states may provide an initial state with sufficient overlap with the desired output for the optimization to be successful. This would combat the trainability issues of such approaches [63,64]. Furthermore, our algorithm could be used to benchmark quantum hardware on strongly-correlated states whenever analytical solutions are known for some expectation values. This can be thought of as a type of application-oriented benchmark [65].
There remain many open questions that would be interesting to explore in future works. A clear next step will be to investigate the optimal strategy with which to compile the unitary gates P k to improve the performance for systems of many magnons [66,67]. Extending our method to open boundary conditions and to other models, such as the Hubbard or Kondo models, is another clear research direction. Additionally, we want to stress that our approach can be used to represent an MPS as a deterministic quantum circuit, which would enable direct preparation of these states on a quantum computer. A Non-unitarity of the R matrices We show here that the R matrices (6) appearing in the Algebraic Bethe Ansatz (ABA) are not unitary in general. We start with where * denotes the complex conjugate. In order for R to be unitary, the following conditions must be satisfied, After some straightforward algebra, the last two conditions are seen to imply with n an integer. However, in general the parameters λ that solve the Bethe equations do not satisfy this requirement.

B The monodromy and transfer matrices
In the main text we have used a matrix that reads in components The YB equation for the matrix R(λ), given in (30), can be expressed in terms of R(λ) as The monodromy matrix T (λ) is a linear map from H a ⊗ H 1 ⊗ · · · ⊗ H N to itself, where H a is an auxiliary space and H j the local quantum spaces (j = 1, . . . , N ). Its definition is [9] T (λ) = R aN (λ)R aN −1 (λ) . . . R a2 (λ)R a1 (λ) . (47) where R aj (u) : H a ⊗ H j → H a ⊗ H j , j = 1, . . . , N . (48) Equation (47) reads in components (49) Notice that the contraction of the matrices R aj in the auxiliary space follows the rule (XY ) a c = X a b Y b c , so that upper and lower indices correspond to row and column indices respectively. For the XXZ model we are considering T (λ) is the 2 × 2 operator matrix Equation (49) can be given the tensor network notation by simply permutating the upper indices of the matrices R, (51) whose graphical representation is The case a = 1 and b = 0 reproduces (4) giving the operator B(λ). Finally, the transfer matrix is defined as the trace of the monodromy matrix (50) in the auxiliary space which yields

C General solution for two magnons
The transformation of the ABA into a quantum circuit can be carried out analytically for the case of two magnons and for any number of sites N .
with R 1 and R 2 parametrized by r 1,2 and s 1,2 respectively. The matrix G 0 (11), which defines the starting point of our algorithm, is given by From now on, we follow the right to left ordering of qubits explained in (25). Notice that this is opposite to the customary left to right order used in Appendix B. Namely . (56) When writing down the entries of an ABC matrix, we will always adopt this convention. Substituting into the recursion relations (16), we find for k ≥ 1. The coefficients c k , d k and e k are assumed to be real and positive using the gauge freedom (20). The new coefficients satisfy while e k is given by The initial conditions for these equations are We obtain the first, 2-qubit unitary of the circuit from equation (18) Equation (17) determines the action of 3-qubit unitaries P k>1 on the state |0 to be The recursion relation (58) defining the coefficients c k is easily solved c k = 1 + |s 2 | 2 + · · · + |s 2 | 2k . (66) Recalling that s 2 is just a phase for one-magnon solutions, we obtain that the first columns of G k and P k |0 for two magnons contain the one-magnon solution (21) and the corresponding P k |0 part of (24).

D The XX model
The XX model is obtained when ∆ = 0. It is a free system and thus implies drastic simplifications with respect to the general XXZ model. The two magnon solution (65) reduces in the XX model to We have also derived the explicit three-magnon solution. The R T basic cell is now , (68) with R 1 , R 2 and R 3 parametrized by t 1,2 , r 1,2 and s 1,2 respectively. We will just sketch the first iterations. The initial matrix G 0 (11) has now dimension 2 × 8. The first step (12) results in a 4 × 8 non-unitary G 1 and the 4 × 4 gate P 1 (64). The second iteration leads to G 2 and P 2 both of dimension 8 × 8. From then on each iteration distill a four qubit unitary, whose action on the rightmost |0 is IBM_Montreal. The first row shows the correlators, while the second row shows their weighted error, defined in Eq. 37. The raw observables were calculated to have a mean weighted error of 0.29, ZNE reduced this error to 0.14 and vnCDR to 0.05.
The construction of R σ is not unique, since has the same effect. Consistently, the YB equation guarantees that R σ = R σ .

G Two-magnon exchange matrices
We construct here the exchange matrices M k implementing the unitary version of the YB equation for The rapidity λ is associated with the matrix R 2 in (54). Following the notation there, we call its entries s 1,2 . The rapidity µ determines the matrix R 2 , with entries r 1,2 . The difference of rapidities λ − µ defines an R matrix whose entries we denote as t 1,2 . The YB equation (33) determines t 1 = s 1 r 1 r 2 1 − r 2 2 + s 2 r 2 , t 2 = s 2 − r 2 r 2 1 − r 2 2 + s 2 r 2

. (81)
From the expressions in Appendix C, we obtain We have generalized (58) to c 2 k (s, r) = c 2 k−1 (s, r) s r * + 1 , subject to the initial condition c 0 (s, r) = 1. Hence c k (s 2 ) = c k (s 2 , s 2 ) coincides with the previously defined coefficients c k , and c k (r 2 ) = c k (r 2 , r 2 ) is assumed real and positive. In deriving the above expression we have used that d k (60) and e k (62) are symmetric under the exchange of s 2 and r 2 . Our main result here is the independence of the exchange matrices M k on the anisotropy.

H Error mitigation details
For detailed analysis of the various mitigation techniques we refer the reader to Refs. [56][57][58]. In particular, the implementations here are very similar to those in Ref. [68]. We use identity insertions [69] to scale the noise by a factor of 3. Therefore, we use noise levels 1, 3 in both ZNE and vnCDR. We note that scaling the noise using identity insertions is not optimal, we expect the implementation of ZNE and vnCDR presented here could be improved upon using pulse stretching. For the training circuits used in vnCDR we replace half of the non-Clifford gates, selecting them randomly and replacing them probabilistically with a Clifford gate as detailed in Refs. [57,58]. For the results taken from IBM_Montreal We use 100 training circuits, and all circuits were run with 32000 shots. Therefore, there is a shot overhead for the implementations of CDR and vnCDR over ZNE by factors of roughly 25 and 50 respectively. For the data obtained from the IBM_Mumbai computer we use 8192 shots with 48 training circuits. Therefore, there is a shot overhead for the implementations of CDR and vnCDR over ZNE by factors of roughly 12 and 25 in this case. These parameters reflect the greatest possible numbers of shots and different circuits that can be run in one job. We note that we repeated the experiments shown in the main text several times on different devices. We present the results where the noise is most stable and well behaved. In all our runs we found vnCDR performs the best on average.
Overall, we find each mitigation method we explore is successful in mitigating the effect of noise. However, even for these small systems it is apparent further techniques should be combined to remove the effect of noise further. We note that our circuit compilation strategies can most likely be improved to reduce circuit depth. Furthermore, a combination of dynamical decoupling [70] and various other error mitigation techniques has been shown to produce accurate observables of interest [71]. It would be interesting to apply a similar approach in the context of ABCs.