\section{Introduction}
\label{sec: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
$\vec{g}$ in a Hamiltonian $H(\vec{g})$ are known as quantum phase transitions
and are fully determined by the ground state of the system
$\ket{\psi_0(\vec{g})}$. The relevant descriptors are the ground state
expectation values of some (local) observables $\hat{O}_i$.
Thus we are interested in characterising the
expectation values of a set of observables $\hat{O}_i$ in this ground state, i.e.
$O_i(\vec{g})=\braket{\psi_0(\vec{g})|\hat{O}_i|\psi_0(\vec{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(\vec{g})$ in that ground state. The
behaviour of these order parameters then gives information about the different
phases the ground state may lie in~\cite{deJongh_1998,Fradkin_1978,Laeuchli_2006,Buchta_2005}.

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)~\cite{Peruzzo2014}. 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
$\braket{\psi(\vec\theta) |H| \psi(\vec\theta)}$ over a family of states
$\ket{\psi(\vec\theta)}$. The hope is that this family is sufficiently
expressive that there is some choice of $\vec\theta$ such that
$\ket{\psi(\vec\theta)}$ 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(\vec g)$ at different parameters
$\vec 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 can be motivated by two intuitions. First, ground states of
quantum systems at a phase transition are characterised by long-range
entanglement~\cite{Sachdev_2011}, and creating complex entangled states will in general
require high-depth VQE circuits. Second, Chen, Gu and Wen~\cite{Chen_2010} 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 \cref{sec:vqe_hardness_in_noninteracting_free_fermion_systems}, we strengthen
this result for translation-invariant non-interacting fermionic models and show
that deep VQE circuits with the Hamiltonian Variational
ansatz~\cite{Wecker2015} 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~\cite{Dreyer_2021}, who analytically find the optimal VQE parameters
and energies for the 1D transverse field Ising model with the Hamiltonian
variational ansatz~\cite{Wecker2015}, 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^{-\lambda 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.~\cite{RocaJerat_2023}
and Jayarama and Svensson~\cite{Jayarama_2022}
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.

\textbf{Our results.}
After briefly reviewing the Variational Quantum Eigensolver method
in~\cref{sec:the-variational-method} and various methods for warm starting
optimisation from previous VQE runs in~\cref{sec:warm_starting_optimisation_with_previous_results}
we present our numerical results in~\cref{sec:results}.
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~\cref{sec:the_tfim} is a spin-1/2 chain, the transverse field Ising model
which has a simple and completely understood phase diagram. Secondly
in~\cref{sec:the_bb_chain}, we considered a spin-1 chain, the
bilinear-biquadratic model~\cite{Laeuchli_2006,Buchta_2005}, with a richer but
still well-understood phase diagram. Finally in~\cref{sec:the_2d_ssh_model}, we
also studied a free fermion model in two dimensions, the 2D SSH model as
studied in~\cite{Obana_2019}. This model is exactly solvable, enabling us to
carry out simulations on a $10\times 10$ lattice (100 qubits).
Additional data for all three models can also be found in~\cref{sec:more_data}.

\begin{table}
  \begin{center}
    \begin{tabular}[c]{l|c|c}
      \multicolumn{1}{c|}{\textbf{}} &
      \multicolumn{1}{c|}{\textbf{Order parameters}} &
      \multicolumn{1}{c}{\textbf{VQE energy derivative}} \\
      \hline
      1D TFIM & \cmark & \cmark \\
      2D TFIM & \cmark & \cmark\\
      1D BBC & \cmark & \cmark / \xmark\\
      2D SSH & \cmark & \xmark \\
    \end{tabular}
  \end{center}
  \caption{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.
  }
  \label{tab:result_overview}
\end{table}

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_{\textnormal{VQE}}(p)$ to detect phase transitions in
\cref{sec:estimating_phase_transitions_via_fits}, 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.

\section{Methods}
\label{sec: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.

\subsection{The variational method}
\label{sec:the-variational-method}

The Variational Quantum Eigensolver (VQE)~\cite{Peruzzo2014, McClean_2016} is a
prominent method for finding ground states $\ket{\psi_0}$ of quantum
Hamiltonians by classically optimising the parameters of a parametric unitary
$U(\vec\theta)$.
The overall goal of VQE is to minimize the objective function
\begin{equation}
  f(\vec\theta) = \braket{\psi(\vec\theta) | H | \psi(\vec\theta)}
       =\braket{\psi_i| U^\dagger(\vec\theta) H U(\vec\theta) |\psi_i},
\end{equation}
where $U(\vec\theta)$ is a family of unitaries parametrised by the classical
parameters $\vec\theta$, $\ket{\psi_i}$ an easily prepared initial state, and $H$
the Hamiltonian whose ground state $\ket{\psi_0}$ we wish to prepare.
By the variational principle $f(\vec\theta)$ is lower bounded by the ground
state energy and can reach that bound only if $U(\vec \theta)$ is sufficiently
expressive and there exists some $\vec\theta^*$ such that
$\ket{\psi(\vec\theta^*)} = \ket{\psi_0}$ for a ground state $\ket{\psi_0}$.
We will call $\ket{\psi(\vec\theta_{\textnormal{min}})}$ the \emph{VQE state},
where $\vec\theta_{\textnormal{min}}$ are the parameters found with minimal
$f(\vec\theta)$. 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 \cref{sec:warm_starting_optimisation_with_previous_results}.

As parametric circuits $U(\vec\theta)$ we use the circuits generated by the
Hamiltonian Variational ansatz (HV)~\cite{Wecker2015}. We split all
Hamiltonians as
\begin{equation}
  H = \sum_{i=1}^k H_i = \sum_{i=1}^k \sum_{j} H_{ij}
  \label{eq:hv_ansatz_hamiltonian_splitting}
\end{equation}
where each $H_i$ consists of mutually commuting, local terms $H_{ij}$. E.g.
for the transverse field Ising model one could choose $H_1 = J \sum_{i=1}^{N-1} Z_i Z_{i+1}$
and $H_2 = h_x \sum_{i=1}^N X_i$ and the local terms would be $H_{1j} = J Z_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
\begin{equation}
  U(\vec\theta) = \prod_{l=1}^p U_{\textnormal{layer}} (\vec\theta_l)
  \quad \textnormal{where} \quad
  U_{\textnormal{layer}}(\vec\theta_l) = \prod_{i=1}^k e^{i \theta_{li} H_i}.
  \label{eq:hv_ansatz_circuit}
\end{equation}
and has $p k$ variational parameters.
If the initial state $\ket{\psi_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.

\subsection{Warm starting optimisation with previous results}
\label{sec: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.~\cite{SmoothParameters_2022} note that the optimal VQE
parameters $\theta_{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
\begin{equation}
  H(\vec g) = \sum_i c_i (\vec g) H_i
\end{equation}
with Hamiltonian parameters $\vec g$ and simple classical functions $c_i (\vec g)$.
This makes the VQE cost function
\begin{equation}
  f(\vec\theta; \vec g) = \braket{\psi(\vec\theta) |H(\vec g)| \psi(\vec\theta)}
                        = \sum_i c_i(\vec g) \braket{\psi(\vec\theta) |H_i| \psi(\vec\theta)}.
\end{equation}
So if we evaluate each $\braket{\psi(\vec\theta) |H_i| \psi(\vec\theta)}$
in a separate run we can easily calculate $f(\vec\theta; \vec g)$ for all values
of $\vec g$ at that $\vec\theta$. We used this property to warm start VQE
at a given $\vec g$ by using the optimal $\vec\theta$ observed in previous runs
at different $\vec g$. Self et al.~\cite{Self_2021} go even further
and propose running multiple instances of VQE at different $\vec 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
\cref{sec:the_effect_of_parameter_extrapolation}. While they worked
well to make $E_{\textnormal{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 $\frac{\Delta E_{\textnormal{VQE}}}{\Delta 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.

\section{Results}
\label{sec: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(\vec\theta) = \braket{\psi(\vec\theta)| H |\psi(\vec\theta) }$
from the full state vector and also exact, analytical gradients $\nabla f(\vec\theta)$
using \texttt{Yao.jl}'s~\cite{YaoFramework2019} automatic differentiation
algorithms. Simulations of the 2D SSH model on up to 100 qubits were carried out using
\texttt{FLOYao.jl}~\cite{FLOYao_2022}, a fermionic linear optics backend for
\texttt{Yao.jl} based on the classical poly-time and space algorithm for
simulating fermionic linear optics circuits~\cite{Terhal_2002,Bravyi_2012}.
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 \cref{sec:warm_starting_optimisation_with_previous_results}.
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
\texttt{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 \cref{sec:more_data} 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.

\subsection{The transverse field Ising model}
\label{sec:the_tfim}
The Hamiltonian of the transverse field Ising model (TFIM) on a graph $G = (V, E)$
is
\begin{equation}
  H_{\textnormal{TFIM}} = J \sum_{(i,j) \in E} Z_i Z_j + h_x \sum_{i \in V} X_i + h_z \sum_{i \in V} Z_i
  \label{eq:tfim_hamiltonian}
\end{equation}
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 \geq 0$ and
fix the bias field strength to be $h_z = \frac{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 $\mathbb{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
$\ket{0 \cdots 0} + \ket{1 \cdots 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.

\begin{figure}
  \begin{center}
    \includestandalone{tfim-sketch}
  \end{center}
  \caption{Sketch of the 1D TFIM Hamiltonian and its splitting into sub-Hamiltonians.
    For the ansatz circuits in
    \cref{eq:hv_ansatz_hamiltonian_splitting} we decompose it into the $ZZ$
    interaction on all {\color{c0} odd edges}, the $ZZ$ interaction on all
    {\color{c1} even edges}, the {\color{c2} transverse field} and the
    {\color{c3} 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 \cref{fig:ssh2d-sketch}.
  }
  \label{fig:tfim-sketch}
\end{figure}

For $|h_x| \ll |J|$ the ground state of the TFIM is the simple product state
$\ket{\psi_{\textnormal{ordered}}} = \ket{1 \cdots 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| \gg |J|$ the ground state is similarly simple:
$\ket{\psi_{\textnormal{disordered}}} = \ket{+ \cdots +}$ and the system is said
to be in the disordered phase.
In 1D, the Kramers-Wannier Duality~\cite{Fradkin_1978} can be used to show that
the phase transition is exactly at $\left|\frac{h_x}{J}\right| = 1$.

As an order parameter to witness the phase transition in the TFIM we use
the $z$-magnetisation
\begin{equation}
  m_z = \frac{1}{N} \sum_{\textnormal{sites } i} Z_i
  \label{eq:z_magnetization}
\end{equation}
which vanishes in the disordered phase and is $\pm 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.

\begin{figure}
  \begin{center}
    \input{figures/tfim_phase_transition_prediction.pgf}
  \end{center}
  \caption{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 $h_x$
    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.
  }
  \label{fig:tfim_phase_transition_prediction}
\end{figure}

\begin{figure}
  \begin{center}
    \input{figures/2dtfim_phase_transition_prediction.pgf}
  \end{center}
  \caption{Predicting the phase transition of the $5 \times 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 $h_x = 3.044$, the phase transition
    point in the thermodynamic limit \cite{deJongh_1998,Bloete_2002}. In the
    lower plot the results for the exact ground state found are shown with red
    dashed lines.
  }
  \label{fig:2dtfim_phase_transition_prediction}
\end{figure}

In~\cref{fig:tfim_phase_transition_prediction,fig:2dtfim_phase_transition_prediction}
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
$z$-magnetisation 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 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 $z$-magnetisation 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~\cite{deJongh_1998,Bloete_2002}.
At first, it might appear surprising that the VQE energy
derivative behaves so similarly before and after the phase transition (remember
that the initial state is the $h_x = 0$ ground state). But it turns out
that evolution with $\sum_i X_i$ and $\sum_i Z_i$
can be combined to apply $R_y(\pi / 2)$ in the first VQE layer, taking
the ground state at $h_x = 0$ to the ground state at $J = 0 = h_z$.
In \cref{sec:more_tfim_data} 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
\cref{fig:tfim_phase_transition_prediction,fig:2dtfim_phase_transition_prediction}
to show the vanishingly small energy derivatives for larger $p$ more clearly.

\subsection{The bilinear-biquadratic chain}
\label{sec:the_bb_chain}

The Hamiltonian of the bilinear-biquadratic chain on $L$ sites is
\begin{equation}
  H_{\textnormal{BBC}} = \cos(\phi) \sum_{i}^{L-1} \vec S_i \cdot \vec S_{i+1}
                       + \sin(\phi) \sum_{i}^{L-1} (\vec S_i \cdot \vec S_{i+1})^2
  \label{eq:bb_hamiltonian}
\end{equation}
where $S_i^{(\alpha)}$ is the spin-1 operator in $\alpha$-direction on site $i$.
The only free parameter in this model is the angle $\phi$.

\begin{figure}
  \begin{center}
    \includestandalone{bbc-sketch}
  \end{center}
  \caption{Sketch of the bilinear-biquadratic chain Hamiltonian. For the ansatz
    circuits we decompose \cref{eq:bb_hamiltonian} into the
    {\color{c0} linear interaction on odd edges},
    the {\color{c1} linear interaction on even edges},
    the {\color{c2} quadratic interaction on odd edges}
    and the {\color{c3} quadratic interaction on even edges}.
  }
  \label{fig:bbc-sketch}
\end{figure}

One checks readily that the terms $\vec S_i \cdot \vec S_{i+1}$ (and hence
also the terms $(\vec S_i \cdot \vec S_{i+1})^2$) commute with the total
spin in $z$ direction
\begin{equation}
  S_{\mathrm{total}}^z = \sum_i S_i^z
\end{equation}
and by $SO(3)$ symmetry the same is true for the analogously defined
$S_{\mathrm{total}}^x$ and $S_{\mathrm{total}}^y$. This makes
$S_{\mathrm{total}}^\alpha$ ($\alpha = x, y, z$) conserved quantities of our
ansatz circuits.

In \cref{fig:bbc_phase_diagram} we show the phase diagram of the
bilinear-biquadratic chain as given in~\cite{Laeuchli_2006,Buchta_2005}. At the
AKLT point $\phi = \arctan (1/3)$ the ground state is exactly expressible as a
matrix product state with bond dimension two, the AKLT state~\cite{AKLT_1987},
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_{\mathrm{total}}^\alpha$ ($\alpha = 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~\cite{Schoen_2005} or using its
symmetries and fusion measurements even in constant depth~\cite{Smith_2022}.

\begin{figure}
  \begin{center}
    \includestandalone{bbc-phase-diagram}
  \end{center}
  \caption{The phase diagram of the bilinear-biquadratic chain. The boundaries
    of the Haldane, trimerised, ferromagnetic can be exactly calculated, but the
    existence of a spin nematic phase in $-3 \pi / 4 < \phi \lesssim -0.67 \pi$ is still
    disputed~\cite{Laeuchli_2006,Buchta_2005}.
  }
  \label{fig:bbc_phase_diagram}
\end{figure}

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 \emph{string order} parameter
\begin{equation}
  O_{\textnormal{string}}^\alpha
  = \lim_{|j-i| \to \infty} S_i^\alpha e^{i\pi (S_{i+1}^\alpha + \cdots + S_{j-1})}   S_j^\alpha
  \label{eq:string_order}
\end{equation}
with $\alpha = x, y, z$. Due to the symmetry of the Hamiltonian, it does not
matter which $\alpha$ one chooses and we used $\alpha = z$ throughout.
For the data shown in \cref{fig:bilbiq_phase_transition_prediction}
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 \emph{staggered dimerisation}
\begin{equation}
  O_{\textnormal{dimer}} = \frac{1}{L-2} \sum_{i=2}^{L-1} (-1)^i (\vec S_{i-1} \cdot \vec S_i - \vec S_{i} \cdot \vec S_{i+1})
  \label{eq:dimerisation}
\end{equation}
which vanishes outside the dimerised phase and is $2$ for a perfectly dimerised
state.
The ferromagnetic phase is indicated simply by the \emph{spin correlation}
\begin{equation}
  O_{\textnormal{spin-corr}} = \frac{1}{L-1} \sum_{i=2}^L  \vec S_{i-1} \cdot \vec S_i
  \label{eq:spin_correlation}
\end{equation}
which again takes the value 2 for a perfectly ferromagnetic state.

\begin{figure}
  \begin{center}
    \input{figures/bilbiq_phase_transition_prediction.pgf}
  \end{center}
  \caption{
    Comparison of the VQE energy derivative (top), string order (second subplot),
    dimerisation (third subplot) and spin correlation (bottom)  as a function of
    $\phi$ 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.
  }
  \label{fig:bilbiq_phase_transition_prediction}
\end{figure}

Our numerical results for the different techniques of detecting phase transitions
for the 12-qutrit BBC are compared in \cref{fig:bilbiq_phase_transition_prediction}.
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 \cref{eq:string_order}, the staggered dimerisation
\cref{eq:dimerisation} and the spin correlation \cref{eq:spin_correlation}
for the same different ansatz depths $p$ and the exact ground state.
Compared to the TFIM and \cref{fig:tfim_phase_transition_prediction} the
situation is less clear for the bilinear-biquadratic chain. The spin-correlation
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.
\cref{fig:bilbiq_vqe_fidelities}) due to the conserved quantities of the ansatz.
The existence of a Haldane phase is correctly detected by the string order
for $p \geq 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 $\phi = - \frac{3 \pi}{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 $\phi$ and simulated system sizes $L$ in \cref{sec:more_bilbiq_data}.

\subsection{The 2D SSH model}
\label{sec:the_2d_ssh_model}

We use the following Hamiltonian for the 2D SSH model
\begin{equation}
  \begin{aligned}
    H_{\textnormal{SSH}} &= -v \kern-1em \sum_{\substack{\textnormal{intercell} \\
                                          \textnormal{edges } \braket{i,j}}} \kern-1em
                     c_i^\dagger c_j + c_j^\dagger c_i
              -w \kern-1em \sum_{\substack{\textnormal{intracell} \\
                                 \textnormal{edges } \braket{i,j}}} \kern-1em
                     c_i^\dagger c_j + c_j^\dagger c_i  \\
          &+ \mu (n_{L} + n_{L^2-L+1} - n_{1} - n_{L^2} )
  \end{aligned}
  \label{eq:2d_ssh_hamiltonian}
\end{equation}
where the $c_i$ are fermionic annihilation operators located on the vertices of
an $L \times L$ square lattice that is subdivided into $2 \times 2$ unit cells
as shown in \cref{fig:ssh2d-sketch}. 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~\cite{Obana_2019}.
Throughout the paper, we fixed the energy scale $v + w = 2$ and the onsite
potential strength to $\mu = \frac{1}{L^2}$. This means the only free parameter
is the hopping strength ratio $\frac{v}{w}$. The addition of $\mu$ ensures that
the ground state in the limit $\frac{v}{w} \to 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}$.

\begin{figure}
  \begin{center}
    \includestandalone{ssh2d-sketch}
  \end{center}
  \caption{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 $-\mu$ and on the top right and bottom
    left corner a positive potential with strength $\mu$ 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 \cref{eq:hv_ansatz_hamiltonian_splitting}.
  }
  \label{fig:ssh2d-sketch}
\end{figure}

In the limit $v \gg 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~\cite{Cade_2019} it
can be prepared in depth $L$ using Givens rotations~\cite{Jiang_2018}. Using
a local encoding~\cite{Klaassen_2021} the same scheme based on Givens rotations
can be used to prepare the trivial ground state in constant depth.
In the limit $v \ll 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 \cref{fig:ssh2d-sketch}.
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 \cref{eq:hv_ansatz_circuit}
for the ansatz circuits. Instead, we use the ansatz circuits already used
in~\cite{Cade_2019} 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 \emph{corner occupation order parameter (COOP)}
\begin{equation}
  O_{\textnormal{COOP}} = e^{i \pi (n_{1} + n_{L} + n_{L^2-L+1} + n_{L^2})}.
  \label{eq:coop}
\end{equation}
In the trivial limit $v \gg w$ the occupation per site in the ground state at half filling
is $\frac{1}{2}$ on all sites, whereas in the topological limit $v \ll w$ the four
corners have occupation $0$ or $1$. This implies that $O_{\textnormal{COOP}} = 0$
in the trivial limit and $O_{\textnormal{COOP}} = 1$ in the topological limit.

\begin{figure}
  \begin{center}
    \input{figures/ssh_phase_transition_prediction.pgf}
  \end{center}
  \caption{
    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
    \times 10$ SSH model. The complete, non-normalized data for the
    VQE energy derivatives is also shown in \cref{fig:ssh_energy_derivatives}.
  }
  \label{fig:ssh_phase_transition_prediction}
\end{figure}

Our numerical findings for the 2D SSH model are summarized in
\cref{fig:ssh_phase_transition_prediction}. Again, we show the VQE energy
derivative in the top subplot and the order parameter---in this case, the corner
occupation order parameter \cref{eq:coop}---in  the bottom. In comparison to
the same data for the TFIM shown in \cref{fig:tfim_phase_transition_prediction}
the situation is quite different 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 \geq 15$ is
surprising---but encouraging in the context of this paper---given the very bad
fidelities shown in \cref{sec:more_ssh_data} between those two states.
The marked difference between the VQE energy derivative in the TFIM (shown
at the top of \cref{fig:tfim_phase_transition_prediction}) 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
\cref{sec:the_tfim}) 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
\cref{sec:vqe_hardness_in_noninteracting_free_fermion_systems} 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.

\section{Discussion and Outlook}
\label{sec:disussion_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
$\left| \frac{E_{\textnormal{VQE}}(p+1) - E_{\textnormal{VQE}}(p)}{E_{\textnormal{VQE}}(p)} \right|$ 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| \sim 20$) to $10^{-3}$ (setting $p=10$ in the 12-qutrit
bilinear-biquadratic chain where $|E_0| \sim 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~\cite{Stanisic_2022,Yu_2022}.

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 $\braket{m_z} \sim \left(\frac{h_x}{J} - 1\right)^{\frac{1}{8}}$ meaning that to
locate the phase transition up to precision $\epsilon$ we need to learn
$\braket{m_z}$ only up to precision $\epsilon^{\frac{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 $\braket{O_{\textnormal{spin-corr.}}}$ as a function of $\phi$
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~\cref{fig:bilbiq_vqe_hardness}
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~\cite{Dreyer_2021}
and~\cref{sec:vqe_hardness_in_noninteracting_free_fermion_systems} 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~\cref{fig:tfim_phase_transition_prediction} 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
~\cite{McClean_2018,Holmes_2021} 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~\cite{Wiersema_2020,Bosse_2022}. 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~\cite{SmoothParameters_2022,Self_2021}. 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~\cite{Dreyer_2021} 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~\cite{Vanhecke_2019},
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~\cref{sec:vqe_hardness_in_noninteracting_free_fermion_systems} we also show
that in the case of translation-invariant free fermion systems 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~\cite{Chen_2010} 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 Hamiltonian.

\begin{acknowledgments}
  The authors would like to thank Filippo Gambetta and other
  members of the Phasecraft team for helpful discussions and feedback on the
  early drafts of this paper.
  This project has received funding from the European Research Council (ERC)
  under the European Union's Horizon 2020 research and innovation programme
  (grant agreement No.\ 817581) and from the EPSRC grant EP/S516090/1.
  All data is  available upon reasonable request.
\end{acknowledgments}