Sketching phase diagrams using low-depth variational quantum algorithms

Mapping out phase diagrams of quantum systems using classical simulations can be challenging or intractable due to the computational resources required to simulate even small quantum systems far away from the thermodynamic limit. We investigate using quantum computers and the Variational Quantum Eigensolver (VQE) for this task. In contrast to the task of preparing the exact ground state using VQE, sketching phase diagrams might require less quantum resources and accuracy, because low fidelity approximations to the ground state may be enough to correctly identify different phases. We used classical numerical simulations of low-depth VQE circuits to compute order parameters for four well-studied spin and fermion models which represent a mix of 1D and 2D, and exactly-solvable and classically hard systems. We find that it is possible to predict the location of phase transitions up to reasonable accuracy using states produced by VQE even when their overlap with the true ground state is small. Further, we introduce a model-agnostic predictor of phase transitions based on the speed with which the VQE energy improves with respect to the circuit depth, and find that in some cases this is also able to predict phase transitions.


I. INTRODUCTION
In a quantum system composed of a large aggregate of particles, a complete description of the full state vector is impractical.A useful approach at the basis of statistical mechanics consists of describing the system in terms of a few properties that describe the collective behaviour of the particles of the system.Abrupt changes in these descriptors as a function of external parameters signal transitions between different phases.At zero temperature, transitions driven by changes of some set of parameters g in a Hamiltonian H(g) are known as quantum phase transitions and are fully determined by the ground state of the system |ψ 0 (g) .The relevant descriptors are the ground state expectation values of some (local) observables Ôi .Thus we are interested in characterising the expectation values of a set of observables Ôi in this ground state, i.e.O i (g) = ψ 0 (g)| Ôi |ψ 0 (g) In computational solid state and condensed matter physics zero-temperature phase diagrams are usually obtained by finding the ground state of the system, using e.g.exact diagonalisation, Monte Carlo methods, or tensor network methods like DMRG or variational methods for PEPS and MERA ansatz states and then measuring the order parameters O i (g) in that ground state.The behaviour of these order parameters then gives information about the different phases the ground state may lie in [1][2][3][4].
Quantum computers could allow the production of families of approximate ground states that are not accessible classically.One method for producing such states, which is particularly well-suited to the noisy intermediate-scale quantum (NISQ) regime, is the Variational Quantum Eigensolver (VQE) [5].VQE is a hybrid quantum-classical algorithm to produce a ground state of a quantum Hamiltonian H via the variational principle.A classical optimiser is used to minimise the expectation value ψ(θ)|H|ψ(θ) over a family of states |ψ(θ) .The hope is that this family is sufficiently expressive that there is some choice of θ such that |ψ(θ) is an approximate ground state of H.
A natural method for applying VQE to approximately computing phase diagrams is then to attempt to find the ground state of H(g) at different parameters g using VQE and then measure the order parameter on the quantum computer in that state.If VQE has produced a high-fidelity approximation to the ground state, the measured value will be close to the true order parameter value.However, even if the approximation produced by VQE is low-fidelity, it is still possible that the measured order parameter will be accurate.This is because order parameters are often local observables that may have the same expectation value in two states, even if the fidelity between those states is low and their global structure is different.
One constraint to this approach is that it requires prior knowledge of the order parameter to be measured, which depends on some knowledge of the physics of the system.We therefore also consider a different kind of order parameter: the rate at which the VQE energy improves as the circuit depth p is increased.Compared to measuring order parameters in the VQE state, this method has the advantage of being model agnostic, as one does not need to understand the physics of the model in question to use the correct order parameters, and the only data required are the VQE energies for different circuit depths throughout the phase diagram.We also test the limitations of this procedure in cases where accessing local information cannot distinguish different phases, like in the case of topological phase transitions.This approach is in spirit similar to the work of Mondaini et al. [6] who use the average sign in quantum Monte-Carlo methods as a model agnostic indicator of different phases.We further motivate our approach by two intuitions.First, ground states of quantum systems at a phase transition are characterised by long-range entanglement [7], and creating complex entangled states will in general require high-depth VQE circuits.Second, Chen, Gu and Wen [8] have shown that symmetry-preserving, local low-depth circuits can map ground states of local Hamiltonians to ground states of other local Hamiltonians, if and only if there is a symmetry-preserving adiabatic path connecting the two Hamiltonians.In appendix D, we strengthen this result for translation-invariant non-interacting fermionic models and show that deep VQE circuits with the Hamiltonian Variational ansatz [9] are required to prepare ground states at or near criticality.We also show that the exact circuit depth for ground states near criticality depends on the gap of the Hamiltonian and hence deeper circuits are required closer to the critical point.
Further support for this approach is provided by recent work of Dreyer, Bejan and Granet [10], who analytically find the optimal VQE parameters and energies for the 1D transverse field Ising model with the Hamiltonian variational ansatz [9], and show that when the initial state for VQE and the target Hamiltonian are in the ordered phase, the VQE energy error (the difference between the true ground state energy and the lowest energy found by VQE) scales with the ansatz depth p as e −λp whereas it scales with p −1 if the target Hamiltonian is in the disordered phase and with p −2 at the critical point.In more recent, related work Roca-Jerat et al. [11] and Jayarama and Svensson [12] also study the complexity of different state preparation methods, including variational algorithms, and find that the complexity depends crucially on whether or not the system passes close to a quantum critical point when preparing the state.Related to our work is also that of Okada et.al. [13] who propose optimising low-depth VQE circuits on a classical computer and then using the quantum computer only to measure non-local order parameters to identify topological phases.
Our results.After briefly reviewing the Variational Quantum Eigensolver method in section II A and various methods for warm starting optimisation from previous VQE runs in section II B we present our numerical results in section III.We compare the approach of measuring order parameters in the VQE state and the approach of analysing the VQE energy derivative to detect and locate phase transitions using VQE using a number of well-understood models.First in section III A is a spin-1/2 chain, the transverse field Ising model which has a simple and completely understood phase diagram.Secondly in section III B, we considered a spin-1 chain, the bilinear-biquadratic model [3,4], with a richer but still well-understood phase diagram.Finally in section III C, we also studied a free fermion model in two dimensions, the 2D SSH model as studied in [14].This model is ex- Comparison of the performance of the different methods used to locate phase transitions across the different models studied by us.A check mark means that the corresponding method signalled phase transitions via an extremum, discontinuity or discontinous first derivative.Because in the case of the bilinear-biquadratic chain the VQE energy derivative detected some, but not all phase transitions we put a check mark and cross on the corresponding entry of the table.Due to finite size effects it is hard to further quantify the success.
actly solvable, enabling us to carry out simulations on a 10 × 10 lattice (100 qubits).Additional data for all three models can also be found in appendix A.
We find that both methods -searching for the ground state using VQE and then measuring order parameters, and analysing the VQE energy derivative -are viable ways to identify phase transitions in the systems we studied.However, neither method is effective for all phase transitions in all systems we studied.For the transverse field Ising model, both methods correctly predict the phase transition and the VQE energy derivative seems to be less prone to finite-size errors than the magnetisation.For the bilinear-biquadratic chain, the VQE energy derivative correctly locates three of the five phase transitions while the order parameters corresponding to the different phases correctly detect those, but have marked finite-size errors in the small system sizes that we simulated.For the 2D SSH model the behaviour is less conclusive and the VQE energy derivative fails to give a clear signal at the phase transition, while an appropriately chosen order parameter measured in the VQE state correctly identifies the phase transition.This indicates that the VQE energy derivative cannot capture topological phase transitions.Besides the VQE energy derivative we also experimented with fitting different functions to the curve E VQE (p) to detect phase transitions in appendix C, but found that this gave no clear advantage.
Crucially, the VQE energy derivative and order parameters were both able to correctly identify and locate phase transitions even when the overlap of the VQE state with the true ground state is very small.This hints that VQE may be able to give interesting information about a physical model, even when the possible circuit depths are not sufficient to prepare the ground state with high fidelity.But neither method works reliably for all models, making the situation analogous to that of trying to detect and locate phase transitions using existing techniques: It requires physical understanding of the model in question and experimentation with different methods to get good results.Nevertheless, our work provides evidence that near-term quantum computers and VQE could be used to understand the phase diagrams of classically intractable models.

II. METHODS
In this section, we give a brief overview of the Variational Quantum Eigensolver and describe our ansatz circuits and the motivation behind them in more detail.We also found that reusing information obtained at different ansatz depths or different points in the phase diagram was crucial to keep the number of VQE runs tractable, and we will describe various methods to reuse this information to warm start optimisation.

A. The variational method
The Variational Quantum Eigensolver (VQE) [5,15] is a prominent method for finding ground states |ψ 0 of quantum Hamiltonians by classically optimising the parameters of a parametric unitary U (θ).The overall goal of VQE is to minimize the objective function where U (θ) is a family of unitaries parametrised by the classical parameters θ, |ψ i an easily prepared initial state, and H the Hamiltonian whose ground state |ψ 0 we wish to prepare.By the variational principle f (θ) is lower bounded by the ground state energy and can reach that bound only if U (θ) is sufficiently expressive and there exists some θ * such that |ψ(θ * ) = |ψ 0 for a ground state |ψ 0 .We will call |ψ(θ min ) the VQE state, where θ min are the parameters found with minimal f (θ).Due to the existence of many local minima, this may not be the global minimum but only the best parameters found with the methods described in section II B.
As parametric circuits U (θ) we use the circuits generated by the Hamiltonian Variational ansatz (HV) [9].We split all Hamiltonians as where each H i consists of mutually commuting, local terms H ij .E.g. for the transverse field Ising model one could choose X i and the local terms would be H 1j = JZ j Z j+1 and H 2j = h x X j .The variational circuit of depth p consists then of p identical layers with k parameters each, takes the form (3) and has pk variational parameters.If the initial state |ψ 0 is the ground state of one of the H i and p is chosen sufficiently large, this variational circuit can represent Trotterised, adiabatic annealing from this H i to the full Hamiltonian H.This means, that for sufficiently large p (and some mild technical constraints on the H i ) this circuit is guaranteed to be able to produce the ground state of H by virtue of the adiabatic theorem of quantum mechanics.However, note that usually VQE is used with relatively small p, where the adiabatic theorem does not hold.

B. Warm starting optimisation with previous results
The cost functions landscapes in VQE usually have many different local minima and finding the global one is hence hard.In all our experiments we want to run VQE for many different ansatz depths p and at many different points in the phase diagram.This makes it possible to warm start the VQE optimisation process by reusing information obtained in previous runs at lower depths or other points in the phase diagram.
Mele et al. [16] note that the optimal VQE parameters θ * li at fixed i often form a smooth function as a function of l.This means one can use good parameters found at low p and warm start optimisation for higher p by using a spline interpolation from the parameters found at low p.However, we found this to be true for some of the models we studied, but not all.For models where it does not hold, one may still warm start the optimisation reusing the optimal parameters found at low p and padding them with zeros to the required length at larger p.
Similarly, one may also want to reuse information obtained at other points in the phase diagram to warm start optimisation.When studying phase diagrams the family of Hamiltonians of interest is of the form with Hamiltonian parameters g and simple classical functions c i (g).This makes the VQE cost function (5) So if we evaluate each ψ(θ)|H i |ψ(θ) in a separate run we can easily calculate f (θ; g) for all values of g at that θ.We used this property to warm start VQE at a given g by using the optimal θ observed in previous runs at different g.Self et al. [17] go even further and propose running multiple instances of VQE at different g in parallel and sharing information between them.But because we use LBFGS instead of Bayesian optimisation as a classical optimiser this information sharing was not possible in our case.
The effect of these two methods to find better initial parameters and reuse already available information is demonstrated in appendix B. While they worked well to make E VQE monotonously decay as a function of ansatz depth p and smoother as a function of the Hamiltonian parameters, we strongly believe that-in particular for large p-they are still not sufficient to find the global minimum.Especially, derived quantities like

∆EVQE ∆p
where still noisy as the VQE energy improves only by very little for large p and small deviations from the global minimum have a large impact here.

III. RESULTS
For all models, we carried out extensive numerical simulations of VQE.For this work the focus was on investigating the optimal case performance of VQE to detect phase transition, so we leave the challenges associated with optimising VQE parameters on real hardware due to noise for future work.Hence we computed exact expectation values f (θ) = ψ(θ)|H|ψ(θ) from the full state vector and also exact, analytical gradients ∇f (θ) using Yao.jl's [18] automatic differentiation algorithms.Simulations of the 2D SSH model on up to 100 qubits were carried out using FLOYao.jl[19], a fermionic linear optics backend for Yao.jl based on the classical poly-time and space algorithm for simulating fermionic linear optics circuits [20,21].For each model, we ran an initial 5 − 10 VQE instances with random parameters initialized uniformly in [0, 1/p].This data was then used to do two more runs with warm starting the optimisation with the previously found parameters, as described in section II B. The first run was warm started with initial parameters extrapolated from better lower depth parameters, if we happened to have found such.The second and final extra run was warm started with initial parameters from a better run at different Hamiltonian parameters, if such a run happened.Unless otherwise specified, all data in this section is from this final run, since it achieved the lowest energy expectation value among all runs.For comparison, we also computed the exact ground state using the Lanczos algorithm implementation in Arpack for the TFIM and the bilinear-biquadratic chain and using full diagonalisation of the single particle Hamiltonian for the 2D SSH model.The overlaps between these exact ground states and the states found by VQE are presented in appendix A for comparison.For the 2D TFIM and the bilinear-biquadratic chain these overlaps mark the phase transitions even more clearly than the order parameters or VQE energy derivatives.But since these overlaps would not be accessible in a real experiment, we did not include them in the main part.For the ansatz circuits in eq. ( 2) we decompose it into the ZZ interaction on all odd edges, the ZZ interaction on all even edges, the transverse field and the bias field.For the 2D TFIM the same decomposition applies, but instead of splitting the interaction part into 2 it is split into 4 parts according to the same edge colouring also used in fig. 7.

A. The transverse field Ising model
The Hamiltonian of the transverse field Ising model (TFIM) on a graph G = (V, E) is where Z i and X i are the spin-1/2 operators in z-and x-direction respectively.Throughout this paper, we set w.l.o.g.J = −1, h x ≥ 0 and fix the bias field strength to be h z = 1 |V | 2 .This means that the only free parameter is the transverse field strength h x (or equivalently the ratio h x /J).The addition of h z breaks the Z 2 -symmetry and ensures that even for finite chains the ground state in the ordered phase is a simple product state-and not the unphysical superposition |0 • • • 0 + |1 • • • 1 -while simultaneously only minimally altering the physics of the system because the gap to the higher energy states scales as |V | −1 .As a simple product state, this initial state can be prepared easily in depth 1.
For |h x | |J| the ground state of the TFIM is the simple product state |ψ ordered = |1 • • • 1 and the TFIM is said to be in the ordered phase.This is also the state we use as an initial state for VQE with the HV ansatz.For |h x | |J| the ground state is similarly simple: |ψ disordered = |+ • • • + and the system is said to be in the disordered phase.In 1D, the Kramers-Wannier Duality [2] can be used to show that the phase transition is exactly at hx J = 1.As an order parameter to witness the phase transition in the TFIM we use the z-magnetisation which vanishes in the disordered phase and is ±1 in the ordered limit, −1 with our choice of a positive bias field strength h z .In the thermodynamic limit, its derivative with respect to the field strength h x is known to diverge at the phase transition.In figs. 2 and 3 we compare our two methods of predicting phase transitions from VQE results for the TFIM.In the 1D case, we see that both the derivative of the zmagnetisation and the VQE energy derivative correctly identify the phase transition and locate it near h x = 1 by having a pronounced peak there.Both methods overestimate the location of the phase transition for small p and converge to lower values with increasing p.However, the argmax of the derivative of the magnetization converges The grey vertical lines are at hx = 3.044, the phase transition point in the thermodynamic limit [1,22].In the lower plot the results for the exact ground state found are shown with red dashed lines.
to a value below h x = 1-most likely due to finite size effects-while the argmin of the VQE energy derivative does not converge to a clear value, most likely due to the derivatives being vanishingly small for large p and hence hard to compute using finite differences from noisy values.The multiple minima for p = 19 in the VQE energy derivatives are most likely spurious and due to the fact that we used finite differences on noisy data to compute the derivatives.
In the 2D case, the VQE energy derivative works well for very low p but is inconclusive for larger p due to the challenges of finding the global minimum of the cost function and the VQE energy derivative being sensitive to this.Unsurprisingly, the derivative of the zmagnetisation works better for estimating the phase transition at higher p.However, due to finite size effects its maximum is not at the point it would be in the thermodynamic limit [1,22].At first, it might appear surprising that the VQE energy derivative behaves so similarly be- FIG. 4. Sketch of the bilinear-biquadratic chain Hamiltonian.
For the ansatz circuits we decompose eq. ( 8) into the linear interaction on odd edges, the linear interaction on even edges, the quadratic interaction on odd edges and the quadratic interaction on even edges.
fore and after the phase transition (remember that the initial state is the h x = 0 ground state).But it turns out that evolution with i X i and i Z i can be combined to apply R y (π/2) in the first VQE layer, taking the ground state at h x = 0 to the ground state at J = 0 = h z .In appendix A 1 we also present the fidelity of the VQE state with the exact ground state for all available data points as well as all data going into the top subplot of figs. 2 and 3 to show the vanishingly small energy derivatives for larger p more clearly.

B. The bilinear-biquadratic chain
The Hamiltonian of the bilinear-biquadratic chain on L sites is where S (α) i is the spin-1 operator in α-direction on site i.The only free parameter in this model is the angle φ.
One checks readily that the terms S i • S i+1 (and hence also the terms (S i • S i+1 ) 2 ) commute with the total spin in z direction and by SO(3) symmetry the same is true for the analogously defined S x total and S y total .This makes S α total (α = x, y, z) conserved quantities of our ansatz circuits.
In fig. 5 we show the phase diagram of the bilinearbiquadratic chain as given in [3,4].At the AKLT point φ = arctan(1/3) the ground state is exactly expressible The boundaries of the Haldane, trimerised, ferromagnetic can be exactly calculated, but the existence of a spin nematic phase in −3π/4 < φ −0.67π is still disputed [3,4].
as a matrix product state with bond dimension two, the AKLT state [23], and is four-fold degenerate due to the two spin-1/2 degrees of freedom at each end of the open boundary AKLT state.We use the AKLT state as the initial state for VQE and to enable our ansatz to produce states with different S α total (α = x, y, z) we include the two spin-1/2 degrees of freedom at either end of the chain in the variational parameters.As an MPS with bond dimension 2 the AKLT state can be prepared in linear depth with the use of one ancilla qubit [24] or using its symmetries and fusion measurements even in constant depth [25].
There are different order parameters to detect the different phases of the bilinear-biquadratic chains.The Haldane phase is heralded by the non-vanishing of the string order parameter with α = x, y, z.Due to the symmetry of the Hamiltonian, it does not matter which α one chooses and we used α = z throughout.For the data shown in fig.6 we used i = 2 and j = L − 1 as the best approximation of "in the bulk" we can do with the small systems we are able to simulate.The dimerisation of the chain is measured by the staggered dimerisation which vanishes outside the dimerised phase and is 2 for a perfectly dimerised state.The ferromagnetic phase is indicated simply by the spin correlation which again takes the value 2 for a perfectly ferromagnetic state.Comparison of the VQE energy derivative (top), string order (second subplot), dimerisation (third subplot) and spin correlation (bottom) as a function of φ for different VQE depths p in the 12 qutrit bilinear-biquadratic chain.The grey vertical lines denote the exactly known phase transitions in the thermodynamic limit.For comparison, we also show the results for the exact ground state with red dashed lines.
Our numerical results for the different techniques of detecting phase transitions for the 12-qutrit BBC are compared in fig.6.In the top subplot, we show the VQE energy derivative, again normalized to have a maximum absolute value of 1 for all different p.In the three lower subplots we show the string order eq.( 10), the staggered dimerisation eq. ( 11) and the spin correlation eq. ( 12) for the same different ansatz depths p and the exact ground state.Compared to the TFIM and fig. 2 the situation is less clear for the bilinear-biquadratic chain.The spincorrelation signals the ferromagnetic phase well, even for low p and even though the true ground state can never be reached using our VQE ansatz (cf.fig.13) due to the conserved quantities of the ansatz.The existence of a Haldane phase is correctly detected by the string order for p ≥ 7 here, similar to the dimerised phase which is detected by the staggered dimerisation.But in both cases, the exact point of the phase transition cannot be detected, partially due to finite size effects as shown by the data for the exact ground state.The VQE energy derivative marks some of the phase transitions with discontinuities but behaves differently depending on the VQE depth.For some p the VQE energy derivative correctly signals the ferromagnetic phase, for others the dimerised or trimerised phase and the AKLT point is clearly marked as the point where the VQE energy stays constant, simply because here the initial state is already the ground state.However, the phase boundaries of the AKLT phase are not clearly marked by the VQE energy derivative.This makes sense because the AKLT phase is a topological phase and as such, it cannot be detected by local observables.The situation around φ = − 3π 4 and whether a phase between the ferromagnetic and the dimerised phase exists is unclear from all our methods.For comparison, we also show the VQE energy derivative for all depths p, parameters φ and simulated system sizes L in appendix A 2.

C. The 2D SSH model
We use the following Hamiltonian for the 2D SSH model where the c i are fermionic annihilation operators located on the vertices of an L × L square lattice that is subdivided into 2 × 2 unit cells as shown in fig. 7. The total number of electrons is L 2 /2, so the system is at half-filling.The second line are positive onsite potentials on the upper right and lower left sites and negative onsite potentials on the upper left and lower right site.Up to the onsite potentials in the second line, this is the same model studied in [14].Throughout the paper, we fixed the energy scale v + w = 2 and the onsite potential strength to µ = 1 L 2 .This means the only free parameter is the hopping strength ratio v w .The addition of µ ensures that the ground state in the limit v w → 0 has definite occupation numbers on all four corners while only minimally altering the physics because the excited modes have a minimal energy scaling with L −1 .
In the limit v w the system is in the trivial phase and at half-filling the ground state is the product of putting two particles into the two lowest energy modes in each unit cell.The circuit depth to prepare this initial state depends on the encoding used to map fermions to qubits.Using the efficient Hamiltonian Variational ansatz from [26] it can be prepared in depth L using Givens rotations [27].Using a local encoding [28] the same scheme based on Givens rotations can be used to prepare the trivial ground state in constant depth.In hopping strength is w.On the top left and bottom right corner we added a negative potential and with strength −µ and on the top right and bottom left corner a positive potential with strength µ to make sure the ground state in the topological limit has occupation number 0 or 1 on the corners.Again, we decompose the Hamiltonian into five sub-Hamiltonians according to the five colours we use for the different terms in this sketch to create the ansatz circuits as in eq. ( 2).
the limit v w the system is in the topological phase and the ground state is the product of putting two particles into the two lowest energy modes in each orange square in fig. 7.In absence of the onsite potential the four corners now host topologically protected zero modes, two of which will be filled.Because the salient feature of the 2D SSH model-the absence or presence of these zero modes-only plays a role at half filling we only consider that case from now on.Because we use the Jordan-Wigner transformation to map the fermionic system to qubits, we cannot directly use the ordering of eq. ( 3) for the ansatz circuits.Instead, we use the ansatz circuits already used in [26] to map the fermionic sites on a square lattice to qubits laid out in a square lattice.
As an order parameter to detect the phase transition in the 2D SSH model, we use the corner occupation order parameter (COOP) In the trivial limit v w the occupation per site in the ground state at half filling is 1  2 on all sites, whereas in the topological limit v w the four corners have occupation 0 or 1.This implies that O COOP = 0 in the trivial limit and O COOP = 1 in the topological limit.
Our numerical findings for the 2D SSH model are summarized in fig.8. Again, we show the VQE energy derivative in the top subplot and the order parameter-in this case, the corner occupation order parameter eq. ( 14)in the bottom.In comparison to the same data for the TFIM shown in fig. 2 the situation is quite differ- ent here: While the corner occupation order parameter correctly predicts and locates the phase transition, albeit only for larger ansatz depths p, the sympathetic eye may see a change in behaviour of the VQE energy derivative for lower p around the phase transition.For large p we mostly see numerical noise and the effects of not being able to find the global minimum in the VQE energy derivative.The good agreement between the corner occupation order parameter in the exact ground state and the VQE states for p ≥ 15 is surprising-but encouraging in the context of this paper-given the very bad fidelities shown in appendix A 3 between those two states.The marked difference between the VQE energy derivative in the TFIM (shown at the top of fig.2) and for the 2D SSH is in part explained by the fact that for the TFIM we know how to take the h x = 0 ground state (i.e. the initial state for VQE) to the J = 0 = h z ground state in one VQE layer (see last paragraph in section III A) whereas for the 2D SSH model, the different Chern numbers of the trivial and the topological phase mean that no translation invariant, low depth, non-interacting circuit exists that can take the ground state of one phase to the ground state of the other phase.See appendix D for a proof of this.This is also consistent with the fact that the VQE energy derivative is a local observable and as such not expected to be able to differentiate between the topological and trivial phase, similar to the way the VQE energy derivative didn't mark the boundaries of the AKLT phase in the bilinear-biquadratic chain clearly.

IV. DISCUSSION AND OUTLOOK
We have demonstrated the usefulness of low-depth VQE to map out phase diagrams.Even when the circuit depth is not sufficient to prepare the true ground state with high fidelity, measuring order parameters or the VQE energy derivative is able to detect phase transitions in the systems simulated by us.Overall, measuring order parameters will require less precision in the estimation of observables than measuring the VQE energy derivative, since one needs the raw expectation value of the order parameters and not the difference between two (possibly very close) expectation values.However, considering the VQE energy has the marked advantage of being model agnostic and may prove useful to gain a rough overview of the phase diagram in question.In contrast to measuring order parameters, the VQE energy derivative also has the advantage that it gives more precise answers for lower depths where the energy decreases faster as a function of p whereas for large p noise makes it harder to accurately measure the VQE energy derivative.For order parameters, on the other hand, it is unclear if lower p is advantageous because the light cones of local observables don't span the whole system and finite size effects are hence suppressed, or if larger p give more accurate results because the VQE state has a larger overlap with the (finite size) ground state.
What remains open are questions of experimental feasibility.A quick back-of-the-envelope calculation shows that the relative accuracy needed to estimate the VQE energy derivative accurately enough to accurately locate the phase transition ranges from 10 −4 (setting p = 10 in the 16-site TFIM where |E 0 | ∼ 20) to 10 −3 (setting p = 10 in the 12-qutrit bilinear-biquadratic chain where |E 0 | ∼ 30).These required relative accuracies were the same for an 8-site or 12-site TFIM and 8-qutrit bilinear-biquadratic chain, hinting that they also remain the same for larger system sizes.However, this is still one order of magnitude away from the relative accuracies of 10 −2 achieved with the help of sophisticated error mitigation schemes in [29,30].This is in stark contrast to the situation where one already knows what phases are expected and is able to measure the corresponding order parameter.Now it is often, but not always, sufficient to distinguish between the order parameter being zero outside of the phase of interest, and having a non-zero value inside the phase of interest.Depending on the critical exponent, the order parameter will also very quickly increase close to the phase transition.E.g. in the case of 1D TFIM we have m z ∼ hx J − 1 1 8 meaning that to locate the phase transition up to precision we need to learn m z only up to precision 1 8 .And the phase transitions to the ferromagnetic phase (certainly from the trimerised phase) in the bilinear-biquadratic chain are of first order, meaning there is a finite discontinuity in O spin-corr.as a function of φ that is easy to detect, even from noisy data.
It remains unclear, how the depth needed to accurately locate the phase transitions will scale with the system size.The data shown in fig.14 and similar analysis for the 1D TFIM and 2D SSH model indicate that to see similar behaviour in the VQE energy derivative the needed circuit depth scales with the system size.On the other hand the analysis done in [10] and appendix D show that in the thermodynamic limit for non-interacting fermionic systems only constant circuit depth is required to prepare ground states away from criticality.It is also interesting to see in fig. 2 that the extremum of the VQE energy derivative and magnetisation derivative start out in a value above the point of the phase transition and then due to finite size effects converge with increasing p to a value below it.This could imply that if only small system sizes are feasible low-depth VQE might locate the phase transition more accurately, because the light-cones of local observables don't span the whole system and hence finite size effects are not as strong.
On real hardware the errors present on NISQ quantum hardware and the sampling noise lead to the barren plateau [31,32] phenomenon which will make optimising the VQE cost functions harder.But our specific choice of ansatz circuits was found to exhibit only mild barren plateaus in the case of the 1D XXZ and 1D TFIM chains or the Heisenberg model on the kagome lattice [33,34].Furthermore, the fact that we run VQE for different circuit depths p and at different points in the phase diagram will help to warm start the VQE optimisation and use initial parameters that are already close to the optimal parameters and hence away from the barren plateaus [16,17].One could even consider simulating VQE on very small systems or with low circuit depths on classical computers and use the optimal parameters found there to warm start the optimisation on the real hardware.
Another open direction is to understand our results analytically, as far as possible.Dreyer, Bejan and Granet [10] derive analytical scaling relations for local observables like m z and m x in the finite depth VQE states of the TFIM and show that scaling collapse techniques, similar to the bond dimension scaling in [35], can be applied to order parameters measured in finite depth VQE states to extract critical exponents.However, it is still open whether these techniques also work for non-integrable models.In appendix D we also show that in the case of translation-invariant free fermion sys-tems deep circuits are needed to connect ground states in different phases and that the circuit depth required to prepare the ground state from a trivial initial state depends on the spectral gap of the Hamiltonian; The smaller the gap is the deeper circuits are required to prepare the ground state with high precision.Again, it is open whether these findings also apply to non-integrable models.The results of Chen, Gu and Wen [8] give a partial answer in that short-time evolution with local unitaries can only connect ground states in the same phase.However, they allow for a bigger family of circuits than we do and make no quantitative statements about how the circuit depth depends on the gap of the target Hamil-tonian.
In figs.9 and 10 we show how the fidelity of the VQE state with the exact ground state increases rapidly as we increase the ansatz depth p far away from the phase transition, but increases much more slowly close to the phase transition.The same message can also be seen in figs.11 and 12 where the energy decreases rapidly for low p far away from the phase transition, in particular for large h x , but then converges soon for h x away from 1 while convergence near the phase transition is slower, resulting in a ridge near h x = 1 for the derivative of the VQE energy.Additionally, we can see how for large p the energy derivative landscape becomes less smooth, due to the inability to resolve very small differences in the VQE energy.

Additional data for the bilinear-biquadratic chain
In fig.13 we show infidelity between the VQE state and the ground space for different ansatz depths p and parameters φ.We can clearly see, that except at the AKLT point where the initial VQE state is already the ground state the overlap between the VQE state and the exact ground space never converges to 1, irrespective of how large we allow p to be.This makes it surprising that the different phases are so clearly visible in fig.14 where we compare the VQE energy derivatives of different system sizes.As expected, the VQE energy converges faster for smaller systems.But the qualitative features and signals between different system sizes are the same, although they are better visible for larger system sizes.For L = 5 the VQE energy converged for p ≥ 7 for all φ and any non-zero value of ∆EVQE/∆p is only due to numerical noise, hence we decided to cap the colour scale at 10 −3 .

Additional data for the 2D SSH model
In fig.15 we show the infidelity between the VQE state and the exact ground space for all p and v/w for the 10 × 10 2D SSH model.For most data points it is exactly 1 and only for a few points at p ≥ 20 do we get a fidelity of 10 −5 .For our simulations of the 6 × 6 2D SSH model the fidelities were not quite that bad, for p ≥ 15 we obtained fidelities of 0.2 for v w away from 1.This is in contrast with fig.16 where we see a clear convergence of the energy with large improvements in the VQE energy, particularly in the topological phase at low p and much lower derivatives at larger p.It should be noted that such very low fidelities are generically expected for    The effect of parameter extrapolation in the bilinear-biquadratic chain.In the top subplot we show the final VQE energy errors of the best run as a function of VQE depth p when using initial parameters randomly distributed in [0, 1/p].In the bottom subplot we show the VQE energy errors when using the parameter extrapolation scheme described in section II B. Each line corresponds to a different value of φ, here showing 13 runs evenly spaced in [0, 2π].The line closest to the AKLT point φ = arctan(1/3) is not shown here, because it is a constant at 10 −8 since here the initial state of VQE is already very close to the target state.large systems due to Anderson's locality ctatstrophe.In these cases the fidelity is a too strong measure of similarity between quantum states, since it bounds the error of all possible observables, and not only that of local observables like we considered throughout this paper.
In this section we show that if |ψ i is the manyparticle ground state of a translation invariant, local, non-interacting fermionic Hamiltonian H (1) and we wish to prepare from this the ground state of another such Hamiltonian H (2) via evolution with said Hamiltonians this is possible in constant time iff the Bloch bundles defined by the ground states of H (1) and H (2) are isomorphic.Furthermore we show that if H (1) is trivial in the sense that its eigenstates are simple Fourier modes the time needed to prepare the ground state of H (2) depends on the band gap ∆E of H (2) .
Let's consider two non-interacting fermionic system in d ≤ 3 spatial dimensions and M ∈ N sites per unit cell (or d ≥ 4 spatial dimensions and M = 1 site per unit cell), described by single-particle, translation-invariant Hamiltonians (D4) Now each block H (i) (k) can be diagonalised via some U (i) (k) with eigenvalues E (i) µ (k) s.t.
and thus (D6) These single-particle eigenstates |k, µ (i) (k) can be used to construct the corresponding Bloch bundles E (i) as sub-bundles of the trivial bundle T d × C M with projection π 1 : T d × C M → T d .The base space of the Bloch bundles E (i) is the first Brillouin zone T d k and the fibers are where 1 ≤ N e ≤ M is the number of filled bands and dimension of the vector bundle and π (i) : E (i) → T d are the bundle projection maps.If h (i)mn x−y decays rapidly as a function of x − y, then H (i) (k) mn will be smooth as a function of k.If, furthermore, H (i) (k) has a gap between the N e -th and the (N e + 1)-th eigenvalue, then the fibers π (i)−1 (k) (interpreted as subspaces of C M ) depend smoothly on k and the E (i) are indeed smooth subbundles of the trivial bundle.
A translation invariant circuit that takes the ground state of H (1) to the ground state of H (2) is now given by a family of unitaries parameterised by k that takes the fiber at k in E (1) to the corresponding fiber in E (2) , i.e. a W (k) s.t.W (k)π (1) −1 (k) = π (2) −1 (k). (D8) One possible choice for W (k) is W (k) = U (2) (k)U (1) −1 (k).However since we only need eq.(D8) and not W (k) |µ (1) (k) = |µ (2) (k) a different choice might yield a smoother W (k). If W (k) smoothly depends on k eq.(D8) is equivalent to saying that the induced map W | E (1) : E (1) → E (2) is a bundle isomorphism.In d ≤ 3 or M = 1 such a bundle isomorphism exists if and only if the Chern numbers of E (1) and E (2)  are equal [39].
We can also use this W (k) to construct the effective Hamiltonian whose time evolution takes the groundstate of H (1) to that of H (2) : In the (k, m) basis define it via Γ(k) = i log W (k) ⇔ W (k) = e −iΓ(k) . (D9) Since basis changing (here the inverse Fourier transformation) and taking logarithms commutes we can directly write down Γ in the (x, m) basis as Γ = x,y,m,n mn decays exponentially as a function of x − y.This implies that Γ is local in real space and hence by a Trotter + Lieb-Robinson argument also found in [8] (approximate) time evolution with it can be implemented in constant time.
The circuits produced by the Hamiltonian Variational ansatz are also generated by the evolution with translation-invariant, non-interacting fermionic Hamiltonians.Hence, if no such W (k) generated by the evolution with a short-ranged effective Hamiltonian Γ exists, it is not possible to prepare the ground state of H (2) from |ψ i with low circuit depth p using VQE with the HV ansatz.
As we have discussed, the existence of a short range effective Hamiltonian Γ depends entirely on the analytic properties of W (k). The analyticity of W (k) is related with the existence of an energy gap in the path connecting H (1) and H (2) .This can be shown to be the case if H (1) is trivial in the sense that its eigenstates are simple Fourier modes (i.e.U (1) (k) = id).Then the Wannier functions of H (2) are given by |w x,µ = y,m dk (2π) d e ik(x−y) W (k) µm |y, m .(D11) Theorem 2.6 in [40] then implies that these Wannier functions decay as y, m|w x,µ ≤ C • e −γ|x−y| where γ grows with the band gap ∆E of H (2) .This implies that in this case also the effective Hamiltonian Γ is more strongly localized if H (2) has a large band gap.By the same Trotter + Lieb-Robinson arguments as above this also implies that deep circuits are required to prepare the groundstate of H (2) that have a small gap if we start with a trivial Slater determinant of Fourier modes |ψ i .

FIG. 1 .
FIG.1.Sketch of the 1D TFIM Hamiltonian and its splitting into sub-Hamiltonians.For the ansatz circuits in eq.(2) we decompose it into the ZZ interaction on all odd edges, the ZZ interaction on all even edges, the transverse field and the bias field.For the 2D TFIM the same decomposition applies, but instead of splitting the interaction part into 2 it is split into 4 parts according to the same edge colouring also used in fig.7.

FIG. 2 .
FIG.2.Predicting the phase transition of the 1D TFIM on 16 sites via the measure of the VQE energy derivative (top) and the derivative of the z-magnetization (bottom).The insets show the argmin (argmax) of the respective quantities with respect to the transverse field strength hx as a function of ansatz depth p.In the lower plot the results for the exact ground state found are shown with red dashed lines.

FIG. 3 .
FIG.3.Predicting the phase transition of the 5 × 5 TFIM via the measure of the VQE energy derivative (top) and the derivative of the z-magnetization (bottom).The grey vertical lines are at hx = 3.044, the phase transition point in the thermodynamic limit[1,22].In the lower plot the results for the exact ground state found are shown with red dashed lines.
FIG. 6.Comparison of the VQE energy derivative (top), string order (second subplot), dimerisation (third subplot) and spin correlation (bottom) as a function of φ for different VQE depths p in the 12 qutrit bilinear-biquadratic chain.The grey vertical lines denote the exactly known phase transitions in the thermodynamic limit.For comparison, we also show the results for the exact ground state with red dashed lines.

FIG. 7 .
FIG.7.Sketch of the 2D SSH model.The intra cell (thick lines) hopping strength is v while the inter cell (thin lines) hopping strength is w.On the top left and bottom right corner we added a negative potential and with strength −µ and on the top right and bottom left corner a positive potential with strength µ to make sure the ground state in the topological limit has occupation number 0 or 1 on the corners.Again, we decompose the Hamiltonian into five sub-Hamiltonians according to the five colours we use for the different terms in this sketch to create the ansatz circuits as in eq.(2).

FIG. 8 .
FIG. 8. Comparison of VQE energy derivative (top) and corner occupation order parameter (bottom) as a function of v/w for different VQE depths p in the 10 × 10 SSH model.The complete, non-normalized data for the VQE energy derivatives is also shown in fig.16.

| 2 FIG. 9 .| 2 FIG. 10 .
FIG.9.Infidelity of the VQE state with the ground state as a function of p and hx for the 1D TFIM with chain length L = 16.

FIG. 16 . 4 φFIG. 17
FIG.16.VQE energy derivative as a function of circuit depth p and v/w for the 10 × 10 2D SSH model.The non-smooth values for p ≥ 20 with multiple minima are most likely due to the inability to find the global ground state for all p and v/w and not due to the complex phase diagram of the 2D SSH model.
FIG. 17.The effect of parameter extrapolation in the bilinear-biquadratic chain.In the top subplot we show the final VQE energy errors of the best run as a function of VQE depth p when using initial parameters randomly distributed in [0, 1/p].In the bottom subplot we show the VQE energy errors when using the parameter extrapolation scheme described in section II B. Each line corresponds to a different value of φ, here showing 13 runs evenly spaced in [0, 2π].The line closest to the AKLT point φ = arctan(1/3) is not shown here, because it is a constant at 10 −8 since here the initial state of VQE is already very close to the target state.

H
D1) where x and y are the unit cells of a square lattice in d dimensions and m and n label the M sites in each unit cell.Using the Fourier transform (i)mn (k) |m n| ⊗ |k k| (D3) where the integral goes over the first Brillouin zone T d and H (i)mn (k) = ∆x e −ik∆x h (i)mn ∆x .

Γ
mn x−y |x, m y, n| = x,y,m,n dk (2π) d e ik(x−y) Γ(k) mn |x, m y, n| = x,y,m,n dk (2π) d e ik(x−y) (i log W (k)) mn |x, m y, n| .(D10) If W (k) is smooth as a function of k, then log W (k) is also smooth and by the Paley-Wiener theorem Γ mn x−y = dk (2π) d e ik(x−y) (i log W (k)) Infidelity of the VQE state with the ground space as a function of p and v/w for the 10 × 10 SSH Model.F = ground states|ψ 0 | ψ0|ψVQE | 2 denotes the squared overlap with the whole, possibly highly degenerate ground space here.Note that compared to fig. 9 the colour scale is on a linear scale here and not on a log scale, because the infidelities are much larger for the 2D SSH model.(bottom pane) VQE energy error as a function of p and v/w for the 10×10 SSH Model.The energy error decreases as p increases, although the fidelity barely improves with p.Note also, that the VQE energy error has a maximum as a function of v/w for p > 10 near the critical point v/w.However, noticing this requires a priori knowledge of E0 which is precisely the quantity that one hopes to find using VQE.