The variational principle states that the energy E obtained by using a certain parameterized wavefunction ψ(θ) for a problem Hamiltonian H is a strict upper bound to the ground state energy E ₀ of H: E = ⟨ ψ θ | H | ψ θ ⟩ ≥ E 0 . (1)

The VQE uses the variational principle, Eq. 1, to compute E ₀ of H on a quantum computer. The VQE algorithm is a hybrid quantum-classical algorithm that utilizes resources from quantum and classical computers in an iterative process. The diagram depicted in Figure 1 shows the link between the quantum processing unit (QPU) and the classical processing unit (CPU). The QPU is responsible for carrying out the computation for a certain quantum circuit that generates the state ψ(θ), depending on a set of parameters, and returns the corresponding bitstrings to the CPU obtained after the measurement. The bitstrings are accumulated and processed by the classical unit, fed to an optimizer that suggests the next set of parameters that will lower the energy in successive iterations.

We analyze the Hamiltonian representing the quantum spin model H = ∑ i > j N J i j x x σ i x ⋅ σ j x + J i j y y σ i y ⋅ σ j y + J i j z z σ i z ⋅ σ j z , (2) where N denotes the number of spins, and σ ^x , σ ^y , and σ ^z are the Pauli matrices. Throughout the rest of the paper, we use units such that ℏ = 1 and J’s are dimensionless. If J i j α α = 1 for all i, j = 1, … , N, and α ∈ {x, y, z}, we call H the isotropic anti-ferromagnetic Heisenberg model Hamiltonian. If all coefficients J i j α α are chosen randomly in the interval ( 0,1 ], then we call H the random Hamiltonian. We use both open and periodic boundary conditions and map each spin to a qubit. For most of our simulations we consider bipartite spin lattices with a single exception of a (frustrated) triangular spin lattice.

In general, the final state of a system acted upon by a parametrized ansatz can be written in the form | ψ f θ 〉 = U θ | ψ 0 〉 , (3) where θ are the variational parameters, U is the ansatz, and |ψ ₀⟩ denotes the initial state. In this paper we demonstrate that the following ansatz is sufficient to yield an accurate approximation to the ground state energy. U θ = ∏ l = N − 1 1 ∏ k = N l + 1 U l k θ l k ∏ l = N − 1 1 ∏ k = N l + 1 U k l θ k l , (4) where U p q θ p q = e − i θ p q σ p y σ q x if p = N or q = N , e − i θ p q σ p y σ q x σ N z otherwise. (5)

We elaborate on how to expand Eq. 4. There is an independent index l and a dependent index k. The index l is monotonically decreased from N − 1 to 1. For each value of l, the dependent index k is decreased from N to l + 1. The index l is not decremented until all the values of k are enumerated. These values of l and k describe the indices of the qubits that the operators in Eq. 5 act upon. In every operator, no more than three qubits are involved. The corresponding parameters θ _kl are independent for each unique combination of l and k. The parameters θ _kl affect the phase of the Nth qubit. The number of unitary operators in Eq. 4 is given by N(N − 1). The ordering of the unitary operators is important. In this paper, the specific combination of the operators in Eq. 5 is termed the XY-ansatz. Any other ordering is prone to produce results which may be different from each other. This is consistent with the findings in recent literature [20, 21].

The motivation to keep the number of operators to a maximum of two or three is inspired from the coupled cluster ansatz which can be powerful enough to express relevant states even in this restrictive form [2]. It is then intuitive to try this approach also for the Heisenberg model. Accordingly, such an ansatz is expressed by | ψ f θ 〉 = e − i A θ | ψ 0 〉 , (6) where A can contain sums of products of Pauli operators. The implementation of an ansatz, e.g. the one given by Eq. 6, is not a simple task in general as it requires factorization of the matrix exponential e ^{−iA( θ )} [22]. Factorization creates a series of products of unitary operations, which results in deeper quantum circuits with a large number of gates. In this work, we do not directly implement the ansatz given by Eq. 6, but seek for other ansätze in a factorized form which do not require further factorization. In effect, we create a quantum circuit from a set of operators instead of a sum of operators. Such an approach allows us to build low-depth quantum circuits. Furthermore, from an experimental perspective, it is difficult to build a quantum computing device in which all the qubits work equally well. Some qubits may perform certain gates more efficiently than others. In order to exploit such devices efficiently and to accommodate for experimental imperfections, we proposed the ansatz where all the parameterized gates are placed on only one qubit. All operations of the parametrised gates can be restricted to this single qubit.

For comparison with the XY-ansatz, we consider two different ansätze. The first one is inspired by quantum chemistry. The unitary coupled cluster ansatz restricted to single and double excitations (UCCSD) is shown to produce results with chemical accuracy [23–26]. We consider U θ = ∏ l = N − 1 1 ∏ k = N l + 1 U k l θ k l , (7) where U k l θ k l = ∏ α ∏ β e − i θ k l β α σ k β σ l α if k = N or l = N , ∏ α ∏ β e − i θ k l β α σ k β σ l α σ N z otherwise, (8) where α, β ∈ {x, y, z}. A combinatorial calculation shows that the number of unitary operators in Eq. 7 is given by 3² N!/2!(N − 2)!. Although the total number of terms scales polynomially rather than exponentially, further reductions are always welcome since the redundant terms often slow down the optimization process. Clearly, the operators in Eq. 4 are a subset of those in Eq. 7. The differences between using these two are highlighted in the results section.

The second ansatz is inspired by the problem Hamiltonian. We consider the ansatz which for the one-dimensional lattice Hamiltonians is given by U θ = U p θ p … U 1 θ 1 , (9) where U _p is given U p θ p = ∏ k = 1 N U k p θ k p , (10) and periodic boundary conditions are used. Each U _kp ( θ _kp ) is given by U k p θ k p = ∏ α e − i θ k α σ k α σ k + 1 α if k = N or k + 1 = N , ∏ α e − i θ k α σ k α σ k + 1 α σ N z otherwise. (11)

For a lattice of size N in one-dimension, the ansatz in Eq. 9 has p × 3N unitary operators. The variational Hamiltonian ansatz [27] is itself inspired from adiabatic evolution. The idea is that a combination of ansatz and initial parameters that mimics the adiabatic evolution can have a lower initial energy to start the variational optimization. It has been used for solving the Hubbard-Fermi model [28, 29] and a modified Haldane-Shastry Hamiltonian [30].

A good choice for the initial state |ψ ₀⟩ often yields better variational results. In the case of antiferromagnets, a good |ψ ₀⟩ for the bipartite lattices is known to be the Néel state where one sublattice is initialised with spins anti-parallel to the other sublattice. The Néel state, for an even number of spins, is in the magnetic sector of zero magnetization, where the ground state for the one-dimensional isotropic anti-ferromagnetic Heisenberg model is located. For the frustrated lattice, half of the lattice spins are initialised anti-parallel to the other half, without regard to their location in the lattice. The qubits representing the spins for one group are initialised as zeros and those in the other group (anti-parallel) as ones.

We use the massively parallel simulator Jülich Universal Quantum Computer Simulator (JUQCS) [31, 32] to perform operations on the state vector. We also use Qiskit [33] for small problem sizes. In an actual quantum device, the state vector itself is not accessible. Instead, the quantum device will produce an ensemble of bitstrings consisting of 0 and 1s only, from which the expectation values of the observables can be derived. This raises two issues. First, since the number of samples or bitstrings can only be finite, it is not always clear if finite sampling can accurately represent the underlying probability distribution. Second, we need a procedure to measure the individual terms of the Hamiltonian. While the first is an open problem, recent works have developed efficient methods for the second when it becomes a problem [34–38]. Fortunately, both these problems do not hinder finding the ground state energy of the Heisenberg model. First, on an actual quantum device, we do not need explicit knowledge of the probability distribution; one can calculate the expectation values directly by sampling. From the samples we can estimate the energy with an accuracy proportional to the square root of the number of samples. Second, unlike in quantum chemistry, the measurement of individual terms for the Heisenberg model is not a problem.

The implementation of the ansatz is best understood through an example. Consider the four-qubit circuit shown in Figure 2, implementing the XY-ansatz for the N = 4 isotropic Heisenberg ring (see Eq. 2). There are 12 parameters in total, and Figure 2 shows the implementation of the first three. After the initial state preparation, gates are applied to construct the unitary operators. The circuit in the rectangular solid box corresponds to the implementation of the exp ( − i θ 43 y x σ z ⋅ σ z ) operator having a variational parameter θ 43 y x . The terms containing σ ^x and σ ^y are implemented by changing to the appropriate basis. The rectangular dashed box highlights the implementation of the third term in the ansatz, given by exp ( − i θ 32 y x σ 3 y ⋅ σ 2 x ⋅ σ 4 z ) . Although the results for the special case of four qubits can be achieved using an even smaller subset of terms in the XY-ansatz operators containing only five parameters, for consistency and completeness, the twelve parameters are used for all the cases. The simulation for this example using Qiskit and JUQCS gives the energy −8.0, which is also the theoretical value obtained by exact diagonalization. The final state also has a 100% overlap with the ground state. The corresponding z- and total-magnetization were both zero, as required.

For the classical optimisation part of the VQE, we use the sequential least squares quadratic programming (SLSQP) [39] optimizer in the SciPy package [40]. SLSQP is a quasi-Newton gradient based algorithm. When using an ideal simulator, the gradients are computed numerically using the cost-effective forward differences formula. The calculation of the energy for given values of the parameters involves the quantum subroutine of VQE. Once the energies are obtained, the calculation of the gradient is done on the classical computer. Calculation of gradient and the next iterate is not a hard problem and classical computers can be used. It is the calculation of the energy for which quantum resources will be required. In our work, the cost of gradient computation for the optimizer is included in the total energy evaluations. For comparison, we also use a gradient-free optimizer called constrained optimization by linear approximation (COBYLA) [40–42].

The results for the one-dimensional isotropic Hamiltonians with periodic boundary conditions are shown in Figure 3A. The coloured lines show the optimization progress, i.e. the lowest energy achieved by successive energy evaluations for different lattice sizes. Initially, the system is in the Néel state. To take advantage of the Néel state, all the optimization parameters are initialised as zeros. As the optimization proceeds, the drop in energy is visible for all lattice sizes shown. The stepped progression is characteristic of the quasi-Newton optimization algorithms, which calculate the gradient before deciding on the step size. Each “step” of the staircase corresponds to a length equal to the number of parameters, since the formula of the forward differences for calculating the gradients requires N (N − 1) + 1 energy evaluations, where N (N − 1) is the number of parameters. Below each optimization curve is a horizontal black line corresponding to the energy of the ground state, obtained numerically by the Lanczos algorithm. According to the variational principle, the calculated energies are upper bounds to the energy of the ground states. Starting from a lattice with 14 spins, the optimizer gave no signal for termination, but due to a time limit of 24 h on the supercomputer [43], the calculations stopped. Due to the different number of compute nodes, circuit depths, and parameters, the total number of energy evaluations differs for each lattice size. The energy optimization process can be restarted from the last values of the parameters. An example is shown for the lattice with 25 spins, as can be seen from the longer curve resulting from the extra energy evaluations. In all the cases shown, the XY-ansatz produced a final energy E _f , which is close to the ground state energy.

The absolute difference between E _f and the corresponding ground state energy E ₀ per spin for each lattice size is shown in Figure 3B. Results up to the lattice size of 6 spins match exactly with the ground state energy [data not shown in (A)], and the differences show up for N ≥ 7. The points cluster in two groups, one corresponding to lattices with an even and another one with an odd number of spins. We observe that the final energies for the lattices with an even number of spins are lower than for the ones with an odd number of spins. Additionally, we note that the ground states of the systems with an odd number of spins are degenerate. Through the variational principle, the ground state energies are mapped to the global minima of an energy landscape when an ansatz can express the ground state. There may be multiple global minima for degenerate cases. Due to this reason the energy landscape will be different from the even spin cases. Therefore, we conjecture that finding the ground state of degenerate cases using VQE is more challenging. This is also observed in the simulations we performed.

Figure 4A shows the energy optimization progress using two different sets of initial parameters. The legend indicates the lattice size. The lines can be separated into two groups, top (dashed) and bottom (solid). The bottom group corresponds to the case where the parameters were initialised as zeros, thus taking advantage of the Néel state and starting from a lower energy value. The top group of lines corresponds to parameters initialised as randomly distributed values in the interval ( 0,2 π ], not taking advantage of the Néel state and starting from a much higher energy, often close to E = 0. For all lattice sizes, initializing the parameters as zeros yields large drops in energy for the first few iterations, and E _f is close to the ground state energy. After the significant drop, the (relative) progress slows down as the energy landscape becomes flatter near a minimum, and a large number of iterations is required to decrease the energy further. On the other hand, random initializations of the parameters do not yield significantly big drops in the energy in any of the cases. A prohibitively large number of energy evaluations seems required to obtain the same accuracy as when parameters are initialised as zeros. From a practical perspective, starting from random parameters does not appear to be very useful for the current problem.

Figure 4B shows the energy optimization progress using the UCCSD inspired ansatz given by Eq. 7. Similar to plot (A), plot (B) shows the trend corresponding to the two different initializations of the parameters. Initializing all parameters as zeros is observed to have a significant initial drop in the energy, contrary to the random parameters. Interestingly, energy optimization progresses for cases with random initial parameters appear to be slower in proportion to the increasing lattice size, i.e. larger N have a slower drop in energy. This effect appears to be consistent with what is termed in the recent literature as the barren-plateau [44–47]. The larger lattice sizes, which require a larger number of parameters, lead to vanishingly small gradients. Given that vanishingly small gradients appear when one approaches a minimum, it is difficult to determine whether the random parameters landed in a local minimum of the energy landscape or at a barren-plateau. One way to find this out would be to restart multiple times with new sets of random parameters, but given that the barren-plateau effect is something that one aims to avoid, which is possible by initializing all parameters as zeros in this case, we skip such an approach.

In both plots (A) and (B) in Figure 4, for parameters initialised as zeros, E _f is very close to the ground state energy. This is understandable as the terms in the XY-ansatz are a subset of the terms given in Eq. 7. However, the difference between the two lies in the number of energy evaluations required to reach the ground state energy. The energies obtained with the XY-ansatz have the same (or comparable) values as the energies obtained with the UCCSD inspired ansatz, but much less energy evaluations with fewer parameters.

Despite the success of initializing the parameters as zeros, it should be noted that such an approach does not necessarily give the lowest possible energy given the ansatz. We performed random initializations with one hundred sets of random parameters for the smaller lattices. The results are shown in Figure 5A. We count and plot the number of cases in which the energy found by starting the optimiser using random initial parameters obtained a final energy equal to or lower than that obtained when starting from the Néel state. We observe that the cases drop sharpy as the lattice size increases. However, finding even a single energy lower than that found when starting from the Néel state shows that the minimum energy reached by initializing all parameters as zeros, although being very close to the ground state energy, is still only a local minimum and usually not a global minimum.

Figure 5B shows the energy optimization progress using the ansatz given by Eq. 9. The lines can again be divided into two groups. The shorter lines on the left correspond to p = 1 and the long lines on the right to p = 5. The lines for p = 5 are longer simply because in this case there are five times more parameters than for p = 1, and so five times more energy evaluations are required to compute the gradient per iteration. For this ansatz, the parameters need to be initialised randomly since initializing them with zeros (for the tested cases with N ≤ 10) results in the optimizer being trapped in a local minimum at the first iteration. It is observed that simply increasing the number of parameters in the circuit is no guarantee for improving the minimum energy. While for the lattice with 19 spins a lower energy is reached with p = 1 than with p = 5, the opposite is the case for the lattice with 20 spins. While a larger parameter space may facilitate a broader spectrum of states, the energy landscape may be the cause of this observation as local minima may surround the global minimum and impede reaching it. Since a random initialization would not be able to take advantage offered by the Néel state, and would thus require a much larger number of energy evaluations to converge to E _f , we restrict our study to only one such random initialization.

Figure 6A shows the energy optimization progress for the random case of the anti-ferromagnetic one-dimensional ring. Using the XY-ansatz, the parameters were initialised as zeros. Except for the ring with 10 spins, none of the optimization processes signalled convergence, and the optimization could be continued to reduce the energy further if required. A quick drop in the initial energy is also observed for the random case, and all the energies are reasonably close to the ground state energy.

Figure 6B shows the energy optimization progress using the gradient-based optimizer SLSQP (solid lines) and the gradient-free optimizer COBYLA (dashed lines). We only show the results for the even lattice sizes. The results for the odd lattice sizes are similar. The gradient-based method gives lower energies at each energy evaluation and finds a much lower minimum faster than the gradient-free method. This result confirms the commonly accepted notion that for noiseless functions, gradient-based methods are superior.

We apply the XY-ansatz without any changes. Results for an isotropic ladder Hamiltonian with open boundary conditions are shown in Figure 7A. The results show the energy optimization progress for ladders of size 3 × 2 up to size 13 × 2. For example, the 8 × 2 ladder, where the optimizer did not converge even after two runs of continued optimization, reported E _f /N = −2.212 as compared to the ground state energy per spin E ₀/N = −2.229.

Figure 7B shows results for three cases with square lattices and a 5 × 6 (frustrated) triangular lattice. The 4 × 4 and 5 × 5 lattices had open boundary conditions and the final energies were E _f /N = −2.218 and E _f /N = −2.298 as compared to E ₀/N = −2.297 and E ₀/N = −2.351, respectively. The optimized (unoptimized) circuit depths were 805 (1398) and 1939 (3531), respectively. For the frustrated lattice, E _f /N = −1.817 as compared to E ₀/N = −1.986. The results for the frustrated lattice show a slightly larger gap in the energy obtained and the ground state. This could either be explained by assuming that the ansatz is less suitable for the case with the frustrated lattice or that the energy obtained when initializing all parameters as zeros corresponds to a local minimum far away from the global minimum.

The case of the square lattice of size 6 × 6 with periodic boundary conditions poses a specific difficulty with the parameter optimization. The problem is that the XY-ansatz for 36 spins has 1260 parameters, and only one iteration can be performed as the evaluation of the gradient is possible only once within 24 h, the time per job on the supercomputer. Since the quasi-Newton methods require multiple iterations without losing the internal variables that systematically improve the convergence per iteration, the optimization is ill-suited for larger problems that cannot undergo multiple iterations in one run. However, this problem can be avoided by saving the internal variables of the optimizer, but this approach is beyond the scope of this work. Despite such a drawback, we proceeded with the first few iterations without saving the internal variables and were able to see a reasonable drop in the initial energy. The circuit depth after optimizing the circuit was reduced from 7458 to 3985. The final energy was E _f /N = −2.631 as compared to E ₀/N = −2.715.

The energies obtained for the one-dimensional lattices can be fitted to a line for extrapolation, as shown in Figure 8A. The slope of the line that gives the energy per spin in the thermodynamic limit, is −1.7783, which can be compared to the exact value −1.7726, known from the Bethe ansatz [48]. The reported value only differs from the exact value by 5.7 × 10⁻³. Note that while the variational theorem guarantees that the energy obtained is a strict upper bound to the ground state, this is not necessarily the case when estimating the value in the thermodynamic limit by means of a fitted extrapolation.

One measure of an ansatz’s ability to scale up beyond what classical computers can simulate is to predict, given the available data, the expected difference between the VQE and exact energies. Such a calculation can be performed by extrapolating the available data. Figure 8B shows data for four different ansätze, namely the ansatz from Eq. 9 with p = 1 (green triangles), p = 5 (red inverted triangles), the ansatz from Eq. 7 (orange squares), and the XY-ansatz (blue dots). Although some of the orange squares are lower than the corresponding blue dots, indicating that a lower energy was obtained using the ansatz from Eq. 7 compared to using the XY-ansatz, the difference is small, and the XY-ansatz is the preferred choice because its number of parameters is significantly smaller. For both the p values tested, the obtained data for the ansatz from Eq. 9 puts it out of competition with the XY-ansatz. Moreover, there is no clear pattern that may help predict the behaviour for the ansatz from Eq. 9 beyond the available data. Linear fitting is performed for the XY-ansatz, and the slope is equal to 2.0829 × 10⁻² with the intercept −8.6202 × 10⁻². Using the given slope, we predict that for a 100-spin ring, the expected energy per spin will be higher than the ground state energy per spin by a value of approximately 0.02.