Exact Ising model simulation on a quantum computer

We present an exact simulation of a one-dimensional transverse Ising spin chain with a quantum computer. We construct an efficient quantum circuit that diagonalizes the Ising Hamiltonian and allows to obtain all eigenstates of the model by just preparing the computational basis states. With an explicit example of that circuit for $n=4$ spins, we compute the expected value of the ground state magnetization, the time evolution simulation and provide a method to also simulate thermal evolution. All circuits are run in IBM and Rigetti quantum devices to test and compare them qualitatively.


I. INTRODUCTION
In recent years the quantum computing has dived fully into the experimental realm. Control of quantum systems has improved so much that quantum computing devices have become a near term reality. These experimental advances rely on some criteria proposed in the 2000's by DiVicenzo [1]: scalable physical system to characterize the qubits, simple fiducial qubit state initialization, long coherence times (longer than the gate implementation times), universal set of quantum gates and qubit-specific measurement capability. Although there are some candidates that can fulfill the first criterion, the field is still in an early stage of development of this technology, where the improvement of qubits control is crucial to accomplish the others.
Private companies have also joined the field. Since 2016, IBM offers cloud based quantum computation platform [2]. Any user can run quantum algorithms on their two five qubits devices, their 16 qubits device, and their 20 qubits device which is available for hubs and partners. It is not the only company that has launched this kind of service: Rigetti Computing also allows the use of its 19qubits device on the cloud [3]. Although both companies are betting for superconducting qubits, their respective device characterization is not the same: basic gate sets and qubits connectivity are some of the differences. As more quantum devices are appearing, it is important to find some methods to test their quality when running sophisticated quantum algorithms.
Several approaches on computer's quantumness have already been tried. The first published article using an IBM device tested the violation of Bell inequalities by more than two qubits (Mermin inequalities) [4], or a recent article tests if the 16-qubit IBM device can be fully entangled by generating graph states [5]. Other works tried to exploit different few qubit experiments, such as error correcting codes and quantum arithmetics [6].
On the other hand, the community has not forgot Feynman's original aim for the proposal of construction of a quantum computer [7]: the simulation of quantum systems. Many classical techniques have been developed in that direction, for instance quantum Monte Carlo methods [8][9][10] or tensor networks algorithms [11]. However, the first suffer from the well-known sign problem and the second are only efficient for slightly entangled systems [12]. In the end, very strongly correlated quantum systems, such as those displaying frustration, will need a quantum computer to be efficiently simulated [13]. There are some works that propose quantum algorithms to construct arbitrary Slater determinants, both in one and two dimensions, to simulate the dynamics of the ground state of fermionic hamiltonians, in particular the Hubbard model [14,15]. Other proposals introduce the concept of compressed quantum computation, i.e. simulation of n-spin chain using log n qubits [16]. This method has been tested in one of the IBM's quantum computers also simulating the magnetization of the Ising model [17]: the main difference respect to the work proposed in that paper is we have access to the whole energy spectrum, which allows us to simulate time and temperature evolution as well.
In this work, we implement a four-qubit experiment that could be interesting both as a proposal for testing and comparing devices quality and for its implications in condensed matter physics. We perform the exact simulation of a spin chain proposed in Ref. [18] with an Isingtype interaction. The Ising model is one of the most famous exactly solvable models, i.e. those models that are integrable. Actually, the steps to find a quantum circuit that diagonalizes the Ising Hamiltonian follow the same strategy than the analytical solution of the model. Therefore, the method can be extended to other integrable models like the Kitaev-honeycomb model, which a circuit has already been proposed [19]. As we are performing an exact simulation, we have access to the whole spectrum and not only to ground state: time evolution and thermal states can be simulated exactly as well. This provides a new approach in quantum simulation if an exact circuit is found for those non trivial models, such as Heisenberg model, which have an ansatz to be solved. In particular, for one-dimensional spin chains, the Bethe ansatz [20] is the most successful method and several proposals exist to simulate and extend it to two-dimensions using tensor network techniques [21]. As the one-dimensional Ising model has analytic solutions for arbitrary number of spins and the circuit proposed in this paper can be efficiently generalized to larger number of qubits, the methods outlined in this work can be used to benchmark a quantum computer by seeing how this compares against known solutions.
The paper is structured as follows. In section II we describe the method proposed in Ref. [18] to construct an efficient circuit that diagonalizes the Ising Hamiltonian: the number of gates scales as n 2 and the circuit depth as n log n. In section III we explain briefly the basic concepts of time evolution in quantum mechanics and give a specific example to be simulated using the circuit derived in the previous section. In section IV, we propose two methods to simulate the expected value of an operator for finite temperature. Section V summarizes the properties of the three devices used for this work, two from IBM and one from Rigetti, and in section VI we present the results of ground state magnetization and the time evolution of | ↑↑↑↑ state and compare the three devices according to them. Finally, the conclusions are exposed in section VII.

II. QUANTUM CIRCUIT FOR THE ISING HAMILTONIAN
Let's consider the existence of a quantum circuit that disentangles a given Hamiltonian and transforms its entangled eigenstates into product states. This circuit will be represented by an unitary transformation U dis where H is the model Hamiltonian and H is a noninteracting Hamiltonian that can be written as H = i i σ z i . This diagonal Hamiltonian contains the energy spectrum i of the original one and its eigenstates correspond to the computational basis states. Then, we will have access to the whole spectrum of the model by just preparing a product state and applying U dis .
For the case of Ising Hamiltonian, the steps to obtain the U dis quantum gate are based on the analytical solution of the model [22,23]: i) Implement the Jordan-Wigner transformation to map the spins into fermionic modes. ii) Perform the Fourier transform to get fermions to momentum space. iii) Perform a Bogoliubov transformation to decouple the modes with opposite momentum. Thus, the construction of the disentangling gate can be done by pieces: In the following subsections, we derive the quantum gates needed to implement the above transformation.

A. Jordan-Wigner transformation
Let's start with the Ising Hamiltonian with transverse field where λ is the field strength. The second term has been added to cancel the periodic boundary term, σ x n σ x 1 , after the Jordan-Wigner transformation in order to solve the system as it was infinite. This modified Hamiltonian will have finite size effects that become negligible as n grows.
The Jordan-Wigner transformation corresponds to transform the spin operators σ into fermionic modes c [24]: where c j and c † j are the fermionic annihilation and creation operators acting on the vacuum |Ω c , c i |Ω c = 0, and following the anticommutation rules {c i , c j } = 0 and {c i , c † j } = δ ij . After this transformation the Hamiltonian reads In terms of the wave function, Notice that the coefficients ψ i1···in do not change. Then it will not be necessary to implement any gates on the quantum register to perform this transformation. However, for now on we should take into account we are dealing with fermionic modes, so any swap between two occupied modes will carry a minus sign. In terms of quantum gates, this is translated into the use of fermionic SWAP gate (fSWAP) each time we exchange two modes: which corresponds with the usual SWAP gate followed or preceded by a controlled-Z gate (see appendix A).

B. Fourier Transform
The next step to solve the Ising model consists on getting the fermionic modes to momentum space using the well-known quantum Fourier transform For n = 2 m for some integer m, this transformation can be implemented with a log-depth circuit and using at most two-body quantum gates. This method is called fast Fourier transform and consists in two parallel Fourier transformations over n/2 sites, the even and the odd sites [25]: To implement such a transformation we need a combination of a two-qubit gate, a 'beam-splitter' F 2 , and onequbit gate, the 'phase-delay' ω k n , which applies the socalled twiddle-factor e 2πik/n : where the fermionic anticommutation relation has been taken into account in the −1 element of the F 2 matrix. All together, the Fourier transform gate becomes The appendix shows the explicit decomposition of this gate.
After the Fourier transformation, the Hamiltonian becomes which it is not diagonal yet as modes with opposite momentum are still coupled.

C. Bogoliubov transformation
The last step will consist on finding a transformation which mixes the two modes according to This will be implemented by a two-qubit gate which acts over qubits that represent opposite momenta. For the case of Ising model, this gate is and its decomposition in basic gates is shown in the appendix. Then, we have finally arrived to the diagonal Hamiltonian: where ω k = λ − cos 2πk n 2 + sin 2 2πk n .

D. n = 4 spin chain
The explicit circuit for a n = 4 chain is shown in Figure  1. First, we prepare the initial state as the ground state for the diagonal Hamiltonian H: The circuit strategy consists in undoing the steps that diagonalize the Ising Hamiltonian. Thus we first undo Bogoliubov transformation by applying (B n k ) † gates, followed by undoing the Fourier transform using the (F n k ) † gates and finally undo the Jordan-Wigner transformation which, fortunately, does not need from any gate as has been explained in the previous section.
For n = 4, the Bogoliubov modes are ±3π/2 and ±π/2, so we need two Bogoliubov gates. Notice that we have removed the (B 0 ) † gate from the circuit of Figure 1; this gate corresponds with the identity for λ > 1 and exchange qubits in the same state for λ < 1, i.e. |00 → −i|11 and |11 → −i|00 . As the initial state for λ < 1 is the |0001 and (B 0 ) † is applied over the last two qubits, it does not affect this state and we can avoid it. If we want to obtain an excited state which eigenstate in the diagonal basis contains |00 or |11 states, then we should only apply bit flip gates over the last two qubits to implement the (B 0 ) † gate.
The circuit shown in Figure 1 also contains fSWAP gates represented with crosses. These will be necessary if even and odd qubits are not physically connected and, as much, they will increase the total number of gates in n 2 . We can eliminate them if the implementation is done in the ibmqx5 device, which allow us to save up to 16 gates of depth according to IBM gate set, but they are indispensable for the implementation in the other IBM device, ibmqx4, as well as in Rigetti's 19-qubit chip.

III. TIME EVOLUTION
Once we have the U dis circuit, we have access to the whole Ising spectrum by only implementing this gate over the computational basis states. This allows us to perform exactly time evolution, where the characterization of all states is needed.
The time evolution of a given state driven by a timeindependent Hamiltonian is described using the time evolution operator U (t) ≡ e −itH : where |ψ 0 is the initial state and i are the energies of the Hamiltonian states |E i . Then, if |ψ 0 is an eigenstate of H there is no change in time (steady state) and therefore the expected value of an observable O will be constant in time. On the contrary, and if [H, O] = 0, the expected value will show an oscillation in time given by We can take advantage from the fact that the eigenstates of the non-interacting Hamiltonian H are the computational basis states and, as we have solved the model, we also know all energies i . Then, it is straightforward to construct the time evolution of a given state |ψ 0 by only expressing it in the computational basis and adding the corresponding factors e −it i . After that, we only need to implement U dis gate over this state to obtain the time evolution driven by the Ising Hamiltonian.
As example, we compute the time evolution of the expected value of magnetization. We take all spins aligned in the positive z direction as initial state, i.e. | ↑↑↑↑ , which in the computational basis is the |0000 state. First, we have to express this state in the H basis, which using U † dis become with φ = arccos(λ/ √ 1 + λ 2 )/2. Then, we apply the time evolution operator to obtain |ψ(t) : Analytically, from which we can obtain the expected value of magnetization, M = 1 2 σ z .

IV. THERMAL SIMULATION
When a quantum system is exposed to a heat bath its density matrix at thermal equilibrium is characterized by thermally distributed populations of its quantum states following a Boltzmann distribution: where β = 1/(k B T ), Z = i e −β i is the partition function and i and |E i are the energies and eigenstates of the Hamiltonian H. The expected value of some operator O for finite temperature is computed as Simulate thermal evolution according to Ising Hamiltonian is, again, straightforward once we have U dis gate, since it consists on preparing the corresponding state in the H basis and apply U dis circuit. In the case of thermal evolution, |E i states are the states of the computational basis, so no further gates are needed to initialize qubits apart from the corresponding combination of X gates.
At that point, we can perform an exact simulation or sampling. In the first case, we run the circuit to obtain the expected value of the observable taking as initial state all states in the computational basis and average them with their corresponding energies. This is done classically once we have the expected values of each state.
On the other hand, we can perform a more realistic simulation by sampling all states according to Boltzmann distribution. First, we need to prepare classically a random generator that returns one of the computational states following the distribution e −β i . Then, we run the circuit many times and compute the expected value of the operator by preparing as initial state the one returned by the generator each time.
The first method demands more runs of the experiment, N × 2 n , needed for the computation of each expected value. No statistical errors come from the averaging part, as it is done classically. For the second method, with only N runs we will obtain a value for the observable with a statistical error of 1/ √ N . For n = 4 the magnetization can be computed analytically and it is shown in Figure 2. At zero temperature, i.e. β → ∞, the system undergoes to phase transition, from paramagnetic to ferromagnetic, at the corresponding critical point of λ = 1. As temperature increases (β decreases), the system have a critical region around λ = 1 until the temperature is high enough to disorder all spins, independently of the transverse field strength (limit β → 0) [26].

A. IBM Quantum Experience
Since 2016, IBM is providing universal quantum computer prototypes based on superconducting transmon qubits which are accessible on the cloud, both interactively in their webpage (Quantum Composer ) [2] or using a software development kit (QISKit). Currently, three quantum chips are available for the general public: two of 5 qubits, ibmqx2 and ibmqx4, and one of 16 qubits, ib-mqx5 [27]. When we performed the experiments ibmqx2 was offline, so we have only used ibmqx4 and ibmqx5.
All backends work with an universal gate set composed by one-qubit rotational and phase gates and a two-qubit gate, the controlled-X or CNOT gate. The differences between the devices, apart from the number of qubits, come from the qubits connectivity or topology and the role that each qubit plays when applied a CNOT gate (control or target). Figure 3 shows the connectivity of the used devices. Each qubit in the 5-qubit device is connected with other two except the central one which is connected with the other four. Qubits in the 16-qubit device are connected with three neighbors in a ladder-type geometry. Both, the one-directionality of CNOT gate and the qubits connectivity, are crucial for the quantum circuit implementation. If the circuit demands an interaction between qubits that are not physically connected, we should implement SWAP gates which will increase our circuit depth and the probability of errors in our final result. Moreover, each time we need to implement a CNOT gate using as a control qubit a physical qubit which is actually a target, we have to invert the CNOT direction using Hadamard gates which, again, will increase the circuit depth and the error probability.
For our propose, ibmqx5 is the best choice for the implementation of the n = 4 circuit. We can use any of the squares and identify upper qubits as 0 and 2 and lower qubits as 1 and 3: according to the circuit of Figure 1, we will not need to use any fSWAP gates. We should only take into account which qubits are control or target to try to reduce the times that we have to invert the CNOT direction.

B. Rigetti Computing: Forest
At the end of 2017, Rigetti Computing launched a 19qubit processor, 'Acorn', that can be used in the cloud through a development environment called Forest [3]. It includes a python toolkit, pyQuil, that allows the users to program, simulate and run quantum algorithms in a similar way as IBM's QISKit. The chip is made of of 20 superconducting transmon qubits but for some technical reasons, qubit 3 is off-line and cannot interact with its neighbors, so it is treated as a 19-qubit device.
Currently, Rigetti's gate set is formed by three one- qubit rotational gates and a two-qubit gate, controlled-Z. This two-qubit gate has the advantage of bi-directionality as the result is the same independently of which is the control qubit. For that reason, the connectivity of the device shown in Figure 4 does not specify the direction of the two-qubit gate. The qubit topology is very different from IBM's devices: some qubits are connected with three neighbors and others with two in a zigzag-type geometry. Then, we can not do without the fSWAP gates, which means that the circuit depth will be greater than the ibmqx5's. On the other hand, it will be comparable with the ib-mqx4, which also needs from these gates. Figure 5 shows the results of the exact simulation of ground state magnetization for the three devices. All points contain a statistical error of 1/ √ N with N = 1024 which comes from the average over all runs to compute the expected value. The other error sources are discussed qualitatively in the following paragraphs.

VI. RESULTS AND DISCUSSION
The best performance come from the ibmqx5 device. This is an expected result as we do not need from fSWAP gates because the qubits connectivity. On the other hand, Rigetti's device performs better than the ibmqx4, even though the number of gates is very similar.
The simulation approaches better to the prediction for low λ. The explanation could come from how affect the experimental error sources to the magnetization. Assuming that two-qubit gates implementation take several hundreds of ns and single qubit gates around one hundred of ns, errors coming from decoherence are expected to be low, as these times are around 50 µs. On the other hand, errors coming from the gate implementation are cumulative and probably the most important error source. It is not negligible neither errors coming from qubits readout, which can induce a bit flip.
The analysis of the results become more clear if we The best simulation comes from ibmqx5 device, which is an expected result since the number of gates used is lesser than with the other devices because of qubits connectivity.
look at the exact ground state wave function: where α = λ − √ 1 + λ 2 and N = 2 √ 2 √ 1 + λα. As λ increases, the amplitude for the states proportional to α goes to zero. That means that any error occurring for λ > 1 is dramatic as it will affect the state with higher probability amplitude, the |1111 . Then, any error in that regime will inevitably cause a decrease in magnetization. On the other hand, errors in some states for λ < 1 can be compensated in average for the other elements with the same probability amplitude.
Similar results are obtained for the time evolution simulation. Figure 6 shows the results for the simulation of the | ↑↑↑↑ state magnetization as it was explained in Section III. Since for the preparation of the initial state it is necessary to implement more gates, we only show the results for the ibmqx5 device, which is the one that can afford this extra circuit depth.
As expected from the previous result, points that represent higher magnetization have more error respect the theoretical values. However, it is remarkable that the relations among the different points for different values of transverse magnetic field are proportionally correct. The oscillations take place for lower values of σ z , have lower amplitudes and are a little bit shifted to the left: even though, they cross each other at the corresponding points and increase and decrease proportionally to the exact result. That is a clear indicator that the error sources in the quantum device are systematic, as the result does not depend on the state preparation.

VII. CONCLUSIONS
In this work we have implemented the exact simulation of a one-dimensional Ising spin chain with transverse field in some quantum computers. To do so, we programmed the algorithm proposed in Ref. [18] in IBM and Rigetti quantum chips. We have simulated the expected value of ground state magnetization as well as the time evolution of the state of all spins aligned. We have also provided two methods to compute thermal evolution of some operator using the same circuit: exact simulation or sampling.
The circuit presented allows to compute all eigenstates of the Ising Hamiltonian by just initializing the qubits in one of the states of the computational basis. It is then a implementation of a Slater determinant with a quantum computer. Since the one-dimensional Ising model is an exactly solvable model, which means that we can compute analytically all the states and energies for any number of spins, and the circuit is efficient, the number of gates scales as n 2 and the circuit depth as n log n, it can represent a method to test quantum computing devices of any size. As has been shown, it is also a hard test since this model is strongly correlated and both the simulation of the phase transition surrounding and time evolution require a high qubits control.
The best performance has been obtained with the ib-mqx5 chip, although the error respect to the theoretical prediction is large in the ferromagnetic phase of the model. A possible reason why this chip shows better results than the others comes from the number of gates used in the quantum circuit, as the qubits connectivity in that device allows us to save all the fSWAP gates. On the other hand, Rigetti's chip performs better than the ibmqx4 chip, even though both implemented circuits have the same gate depth.
The ferromagnetic phase is difficult to simulate due to the fact that any error that can induce a qubit bit flip will produce a decrease in magnetization, as can be traced out from the ground state wave function of Eq. (26). However, and taking into account this fact, the time evolution simulation is reasonably good since the expected oscillations for different transverse magnetic field strengths are shifted to the left and have lower amplitude and magnetization, but are also proportional each other as are the theoretical values.
As a final remark, this circuit is also interesting from a point of view of condensed matter physics as specific methods to simulate exactly time and thermal evolution are provided. This can open the possibility of simulate other interesting models: integrable, like Kitaev Honeycomb model [19], or with an ansatz, like Heisenberg model.

Notes
The program used for this work for the IBM quantum devices was awarded with the IBM "Teach Me QISKit" award [28].
Recently, Rigetti computing changed the quantum device to one of 8 qubits. The results shown in this work correspond to the previous device of 19 qubits. the qubits connectivity; this can be easily done using the identity ( In addition, controlled-Z gate could be implemented using two Hadamard gates and a CNOT, as it is also shown in Figure 7. From Rigetti's implementation point of view, controlled-Z gates are part of their basic gate set and, although CNOT gate and H are not, they are included in pyquil language, so the quantum programmer should not care about decompose them in terms of the other gates (except to keep in mind the circuit depth will increase with the use of non-basic gates).

B. Fourier transform
The building blocks of Fourier transform gate are represented by the quantum gate of Eq.(11), which decomposition is shown in Figure 8. Its implementation requires the controlled-Hadamard gate, decomposed in Figure 9. Phase gate Ph is included in Rigetti's pyquil language and it is equivalent to IBM's U 1 gate (see Eq. (24)). In particular, Ph(π/2) and Ph(π/4) correspond to S and T gates respectively.

C. Bogoliubov transformation
Bogoliubov transformation is implemented using B n k gates written in Eq. (14). The explicit decomposition is shown in Figure 10, where the controlled-RX gate (shown in Figure 11 has been decomposed using the methods of Ref. [29]). Rotational gates are part of the basic Rigetti gate set and are equivalent to IBM's gates R X ≡ U 3 (φ = 0, λ = π), R Y ≡ U 3 (φ = λ = 0) and Figure 10: Bogoliubov gate decomposition. Controlled-RX gate needed is shown in Figure 11.

D. Initial state preparation
The ground state of the n = 4 Ising model in the diagonal basis is |0000 for λ > 1 and |0001 for λ < 1. Qubits are always initialized in the |0 state both in IBM and Rigetti devices. Then, to compute the ground state, we only need to perform a bit-flip gate (Pauli-X or X gate), on fourth qubit to initialize the circuit for λ < 1.
The example given for time evolution simulation requires from the preparation as the initial state the one shown in Eq. (20). This can be done by applying a R Y (φ) gate on the first qubit to introduce the φ angle, followed by a phase gate to introduce the evolution phase e 4it √ 1+λ 2 and a CNOT gate between first and second qubits.