# Energy Landscapes for the Quantum Approximate Optimisation Algorithm

Choy Boy  David J. Wales  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom

###### Abstract

Variational quantum algorithms (VQAs) have demonstrated considerable potential
in solving NP-hard combinatorial problems in the contemporary near
intermediate-scale quantum (NISQ) era. The quantum approximate optimisation
algorithm (QAOA) is one such algorithm, used in solving the maximum cut (Max-
Cut) problem for a given graph by successive implementation of $L$ quantum
circuit layers within a corresponding Trotterised ansatz. The challenge of
exploring the cost function of VQAs arising from an exponential proliferation
of local minima with increasing circuit depth has been well-documented.
However, fewer studies have investigated the impact of circuit depth on QAOA
performance in finding the correct Max-Cut solution. Here, we employ basin-
hopping global optimisation methods to navigate the energy landscapes for QAOA
ansätze for various graphs, and analyse QAOA performance in finding the
correct Max-Cut solution. The structure of the solution space is also
investigated using discrete path sampling to build databases of local minima
and the transition states that connect them, providing insightful
visualisations using disconnectivity graphs. We find that the corresponding
landscapes generally have a single funnel organisation, which makes it
relatively straightforward to locate low-lying minima with good Max-Cut
solution probabilities. In some cases below the adiabatic limit the second
lowest local minimum may even yield a higher solution probability than the
global minimum. This important observation has motivated us to develop broader
metrics in evaluating QAOA performance, based on collections of minima
obtained from basin-hopping global optimisation. Hence we establish
expectation thresholds in elucidating useful solution probabilities from local
minima, an approach that may provide significant gains in elucidating
reasonable solution probabilities from local minima.

††preprint: AIP/123-QED

##  I Introduction

The initial setback of implementing practical quantum algorithms utilising the
quantum phase estimation (QPE) architecture onto current-day near
intermediate-scale quantum (NISQ) devices [1, 2, 3], which typically possess
short decoherence times [4] and significant quantum noise [5], has prompted
the rapid development of variational quantum algorithms (VQAs) with shorter
quantum circuit depths [6, 7, 8]. VQAs typically operate within a hybrid
classical-quantum optimisation framework [9], where an initial quantum state
is evolved by a parameterised circuit ansatz on a quantum device. After the
final evolved wavefunction is measured, a classical optimiser evaluates the
cost function from the measurement and subsequently suggests new parameters
that are fed back into the parameterised quantum circuit. This interface
between classical and quantum computers iterates until a suitable convergence
criterion is attained. VQAs have proved to be surprisingly robust in tackling
various sources of noise attributed to NISQ devices, such as decoherence [10]
and depolarisation [11]. These properties are ascribed to the innate
variational nature acting as a parametric lever, which can be flexibly
adjusted even under noisy environments [12]. Coupled with the recent advances
in error-mitigation methods [13, 14, 15], it is anticipated that VQAs will
enable the realisation of practical quantum advantage before the advent of
fault-tolerant devices equipped with error-correction protocols [16].

The quantum approximate optimisation algorithm (QAOA) is a VQA that
monotonically improves in performance as the number of quantum circuit layers
in the ansatz increases [17]. Variational quantum eigensolver (VQE) algorithms
typically utilise the global minimisation of the cost landscape directly as
the solution of interest [18]. In contrast, QAOA uses the resulting final
wavefunction to obtain approximate solutions based on the states in the
computational basis with the highest frequency when measurement after circuit
evolution is carried out. Hence QAOA is a particularly attractive algorithm
for solving combinatorial optimisation problems, such as Max-Cut [19], with
promising applications in portfolio optimisation [20] and chemistry [21, 22].

The expansion of VQAs has also motivated the study of their potential
shortcomings in solving practical large-scale problems from a software
perspective. In particular, for various VQA ansätze there is an exponential
growth in the barren plateau problem as the number of qubits and circuit
layers required to encode a given problem increases [23, 24, 25]. Recently,
more insight has also been gained into the challenges of exploring the cost
landscapes of VQAs that arise from a proliferation in the number of local
minima and other stationary points as the complexity of the problem increases
[26, 27]. However, further analysis of the organisation of the global cost
landscapes of VQAs, and how this structure impacts the quality of the
solutions obtained, is needed, especially for QAOA [28, 29].

Here, we seek to address these gaps in understanding using the well-
established theory and associated computational methodology of molecular
energy landscapes [30, 31]. We characterise the cost optimisation landscapes
of QAOA for various weighted and unweighted graphs in solving the Max-Cut
problem using basin-hopping global optimisation methods [32, 33, 34] to locate
global minima, and discrete path sampling [35, 36] to create connected
databases of minima and the transition states that connect them. Recently,
energy landscape techniques have demonstrated considerable utility for quantum
computing in the analysis of hardware-efficient ansätze for the VQE algorithm
[37], and optimisation of electronic wavefunctions in a combined discrete
space of operators and continuous parameter amplitudes [38]. Our new results
for QAOA show that the solution landscapes below the adiabatic limit generally
possess single-funnelled structures associated with self-organising systems
where locating the global minimum is relatively straightforward [39, 30, 40].
Furthermore, we find that local minima sufficiently close in energy to the
global minimum may also exhibit good solution probabilities for the Max-Cut
problem. In some instances, the second lowest minimum has a higher solution
probability than the global minimum, highlighting the importance of studying
the VQA solution landscape globally. This observation leads us to introduce
metrics that take into account the distribution of minima in evaluating the
performance and robustness of QAOA. We also utilise the convex hull of the
solution space in estimating expectation cutoffs for the location of local
minima with reasonable solution probabilities. We hope that these techniques
can advance the feasibility of implementing QAOA for problems with numerous
local minima in noisy environments.

##  II Methodology

Given an undirected graph $G=(V,E)$, with weights $w_{ij}$ assigned to edges
$e_{ij}\in{E}$ for connected vertices $i,j\in{V}$, the Max-Cut problem seeks
to partition $V$ into two distinct sets such that the sum of weights between
the two sets is maximised. If $w_{ij}=1$ for all $e_{ij}$, then $G$ is said to
be an unweighted graph; otherwise $G$ is a weighted graph. It follows that the
Max-Cut problem can be mapped to a two-spin Ising-type cost Hamiltonian
$\hat{H}_{C}$ corresponding to $N$ implementable qubits:

$\hat{H}_{C}=\frac{1}{2}{\sum}_{e_{ij}\in{E}}{w_{ij}}(Z_{i}\otimes Z_{j}),$
(1)

where the states $\ket{s}=\\{\ket{\alpha},\ket{\beta}\\}^{\otimes{N}}$ encode
the desired solution strings to the Max-Cut problem, with
$\ket{\alpha}=\ket{0}$ if and only if $\ket{\beta}=\ket{1}$, and vice versa.
Thus, the aim of QAOA is to approximate the ground-state energy or the lowest
eigenvalue of $\hat{H}_{C}$ via a suitable ansatz with unitary operator
$\hat{U}(\bm{\theta})$ to evolve an initial state $\ket{\psi_{0}}$, and
subsequently use the final evolved state
$\ket{\Psi(\bm{\theta})}=\hat{U}(\bm{\theta})\ket{\psi_{0}}$ to approximate
$\ket{s}$. This objective can be achieved on a quantum device by performing a
certain number of shots per experiment and measuring all qubits in the
computational basis, taking the state possessing the greatest number of shots
to best approximate $\ket{s}$ for that experiment. We seek to simulate this
procedure classically by considering the probability of measuring the state
$\ket{s}$ in the computational basis, $p(\ket{s})$:

$p(\ket{s})=|\braket{s}{\Psi(\bm{\theta})}|^{2}.$ (2)

The objective function to be minimised by the classical computer is the
expectation of $\hat{H}_{C}$, $\langle\hat{H}_{C}\rangle$:

$\langle\hat{H}_{C}\rangle={E({\boldsymbol{\theta}})}={\bra{\psi_{0}}\hat{U}^{{\dagger}}(\bm{\theta})\hat{H}_{C}\hat{U}(\bm{\theta})\ket{\psi_{0}}}.$
(3)

The QAOA ansatz with parameters $\bm{\theta}=\\{\bm{\gamma},\bm{\delta}\\}$
can be assembled as a Trotterised variational schedule, comprising a cost
circuit layer with unitary operator $\hat{U}_{C}(\gamma)$, followed by a mixer
circuit layer with unitary operator $\hat{U}_{M}(\delta)$ up to a circuit
depth $L$:

$\ket{\Psi(\bm{\gamma},\bm{\delta})}={\prod_{l=1}^{L}}\hat{U}_{M}(\delta_{l})\hat{U}_{C}(\gamma_{l})\ket{\psi_{0}},$
(4)

where $\ket{\psi_{0}}=\ket{+}^{\otimes{N}}$ is the state encoding for all
possible partitions of $V$ with equal probability. The cost layer
encapsulating $\hat{H}_{C}$ can be compiled as a sequence of two-qubit
parameterised $R_{zz}$ quantum gates for qubits $q_{i}$ and $q_{j}$, with
$\gamma$ scaled based on the weights of $e_{ij}$:

$\begin{split}\hat{U}_{C}(\gamma)&=e^{-i\gamma\hat{H}_{C}}\\\
&={\prod}_{e_{ij}\in{E}}R_{zz}(-w_{ij}\gamma).\end{split}$ (5)

The mixer layer performs a time-evolution of the mixer Hamiltonian
$\hat{H}_{M}=-\sum_{i=1}^{N}{X_{i}}$, which anticommutes with $\hat{H}_{C}$
and is an eigenfunction of $\ket{\psi_{0}}$. The mixer layer can be realised
as a parallelisation of single-qubit parameterised $R_{x}$ quantum gates:

$\begin{split}\hat{U}_{M}(\delta)&=e^{-i\delta\hat{H}_{M}}\\\
&={\bigotimes}_{i=1}^{N}{R_{x}}(2\delta).\end{split}$ (6)

It has been shown that QAOA conforms to the adiabatic theorem, i.e. for
$L\rightarrow\infty$ the final evolved state $\ket{\Psi(\bm{\theta})}$
converges exactly to the ground state of $\hat{H}_{C}$, and thus gives the
optimal $p(\ket{s})$ [17]. In practice, such an implementation is unfeasible
in the NISQ regime, hence we are interested in considering $L_{ad}$ for a
given system, defined as the minimum number of circuit layers required to
reach the adiabatic limit, assuming that $L_{ad}$ can be attained. As we will
demonstrate in our analysis of the energy landscapes of QAOA, it is also
important to distinguish $L_{ad}$ from $L_{min}$, where $L_{ad}\geq L_{min}$.
Here, $L_{min}$ is the minimum number of layers needed to achieve the maximum
$p(\ket{s})$ in the corresponding global minimum. Local minima with lower
$p(\ket{s})$ may be present due to underparametrisation of the circuit ansatz,
and hence a less thorough exploration of states in the Hilbert space may be
sufficient to obtain a useful solution. We hypothesise that the exponential
increase in the number of local minima is attributable to circuit ansätze with
$1\leq{L}\leq{L_{min}}$ layers. The behaviour of local minima may vary for
$L_{min}<{L}\leq{L_{ad}}$ layers if $L_{min}<{L_{ad}}$, and we observe for
various graphs that the number of local minima may increase first before
decreasing to the adiabatic limit, or instead decrease monotonically.

For each graph considered and $L$, we generate an initial set of minima via
basin-hopping global optimisation [32, 33, 34] using the GMIN program [41].
The analytic gradients of the parameterised rotation gates were calculated via
the parameter-shift rule [42] (see Appendix A for more details):

$\frac{\partial{E({\boldsymbol{\theta}})}}{\partial\theta_{i}}=\frac{1}{2}\left[{E\left({\boldsymbol{\theta}}+\frac{\pi}{2}\bm{e}_{i}\right)-E\left({\boldsymbol{\theta}}-\frac{\pi}{2}\bm{e}_{i}\right)}\right].$
(7)

Local minimisation for the basin-hopping steps employed a limited-memory
Broyden [43], Fletcher [44], Goldfarb [45], Shanno [46] (L-BFGS) procedure
[47, 48] equipped with the Metropolis criterion for accepting/rejecting steps
[49]. The resulting minima were then employed as the starting points for
construction of a kinetic transition network [50, 51, 52]. Discrete path
sampling [35, 36] (DPS) was used via connection attempts for selected minima
pairs with final states $\ket{\Psi(\bm{\theta}_{\mu})}$ and
$\ket{\Psi(\bm{\theta}_{\nu})}$. The doubly-nudged [53, 54] elastic band [55,
56, 57, 58] approach was used to locate candidates for accurate transition
state refinement by hybrid eigenvector-following [59, 60, 61]. The missing
connection algorithm [62] was used to select pairs of minima to fill in the
gaps in incomplete pathways via Dijkstra’s shortest path algorithm [63]
combined with a distance metric based on the state overlap between local
minima, $S_{\mu\nu}$:

$S_{\mu\nu}=1-|\braket{{\Psi(\bm{\theta}_{\mu})}}{{\Psi(\bm{\theta}_{\nu})}}|.$
(8)

Any new minima are added to the database along with the transition states and
connection information. The resulting cost landscapes can be visualised using
disconnectivity graphs, where the local minima are segregated into disjoint
sets for regular thresholds in the energy [64, 39]. In these graphs, the
bottom of a branch corresponds to the energy of a local minimum, and the
branches are joined when the corresponding minima can interconvert via a
pathway below the given threshold. Visualisation of the energy landscape can
be further enhanced by colouring each minimum with the corresponding
probability of finding the correct Max-Cut solution; we find this construction
especially useful in comparing the solution landscapes as $L$ varies.

Figure 1:  (a) (left) Weighted graph $G_{1}$ encoded with four qubits. (b)
(right) Corresponding QAOA ansatz for $G_{1}$. Figure 2:  (a) (left) Contour
plot of $\langle\hat{H}_{C}\rangle$ for graph $G_{1}$ and its QAOA ansatz with
$L=1$ against the parameters $\gamma$ and $\delta$. The solid black line
depicts the minimum of each vertical slice of the contour plot, with the
global minimum (GM, in green), local minimum (LM, in red) and the transition
state connecting them (TS, in black) situated on the pathway that is also
plotted. (b) (right) Corresponding disconnectivity graph of the contour plot
in Fig. 2a, with both the GM and LM coloured based on their respective
probabilities of obtaining the correct Max-Cut solution of
$\ket{\alpha\beta\alpha\beta}$.

As an example, consider the weighted graph $G_{1}$ with four vertices (Fig.
1a), where the Max-Cut problem can be encoded as a four-qubit QAOA ansatz with
varying $L$ (Fig. 1b). We note that although the ansatz cost circuit layer can
be compiled in numerous ways, the arrangement in Fig. 1b is ideal in reducing
the overall circuit depth. This ordering does not require swap gates to
permute non-neighbouring qubits corresponding to their respective edges, which
is an important consideration when transpiling QAOA onto planar
superconducting processors [65]. The cost landscape of the $L=1$ circuit
ansatz features a global and a local minimum connected by a transition state
(Fig. 2a), and the corresponding disconnectivity graph is shown in Fig. 2b,
where the branches are coloured with the probabilities of finding the state
$\ket{\alpha\beta\alpha\beta}$, which corresponds to the Max-Cut solution of
$G_{1}$.

##  III Results

###  III.1 Complete unweighted graphs

We begin by examining the QAOA circuit ansätze for the complete unweighted
graph series $K_{N}$ from $N=3$ to $N=8$ and from $L=1$ to $L=3$, where each
vertex is connected to every other vertex with unit weight. It follows that
the Max-Cut solution of $K_{N}$ is the set of all possible tensor product
permutations of $\lfloor N/2\rfloor$ number of $\ket{\alpha}$ states and
$\lceil N/2\rceil$ number of $\ket{\beta}$ states: thus the total number of
Max-Cut solutions for odd $N$ is $2N!/\\{\lfloor N/2\rfloor!\lceil
N/2\rceil!\\}$, and for even $N$ $N!/\\{2(N/2)!\\}$ solutions.

Table 1:  Number of minima $M$ (top value) and the highest correct Max-Cut probability (HCMP) (bottom value) for graphs $K_{3}$ to $K_{8}$ with varying $L$ obtained from basin-hopping global optimisation. For the case of $K_{8}$ and $L=3$, the HCMP is not equal to the maximum value of 1, and the number of decimal places used for all HCMPs is chosen to be the same as this case for ease of comparison. Graph | $L=1$ | $L=2$ | $L=3$
---|---|---|---
$K_{3}$ | 1 | 1 | 1
1.000000 | 1.000000 | 1.000000
$K_{4}$ | 1 | 1 | 1
0.739106 | 1.000000 | 1.000000
$K_{5}$ | 1 | 4 | 1
0.975990 | 1.000000 | 1.000000
$K_{6}$ | 1 | 23 | 324
0.671340 | 0.994239 | 1.000000
$K_{7}$ | 1 | 37 | 598
0.951350 | 0.999619 | 1.000000
$K_{8}$ | 1 | 46 | 3418
0.629727 | 0.991483 | 0.999997

Figure 3: Disconnectivity graphs of $K_{6}$, $K_{7}$ and $K_{8}$ for $L=2$.

Figure 4: Disconnectivity graphs of $K_{6}$, $K_{7}$ and $K_{8}$ for $L=3$.

Table I summarises the number of minima $M$ and the highest correct Max-Cut
probability (HCMP) from the collections of minima obtained via basin-hopping
global optimisation for $L=1$ to $L=3$. The complete graphs with odd $N$
generally possess higher HCMPs than their even counterparts, mainly due to
their greater number of accepted Max-Cut solutions that contribute to their
corresponding probabilities. We also find that although for $K_{6}$ to $K_{8}$
the expected exponential increase in $M$ is observed as $L$ increases, for
$K_{5}$ there was a decrease in the number of minima from $L=2$ to $L=3$,
leading to a simplification in the energy landscape from $L_{min}=2$ to
$L_{ad}=3$. Looking at the disconnectivity graphs of $K_{6}$ to $K_{8}$ for
$L=2$ (Fig. 3) and $L=3$ (Fig. 4), we observe that the majority of the local
minima generally possess very high correct Max-Cut probabilities, particularly
those closer to the global minimum, especially as $L$ and $M$ increase. The
well-funnelled organisation of the landscape also becomes more apparent as $L$
increases, and this structure is expected to outweigh the challenges
associated with solving the Max-Cut problem for higher $N$, particularly for
$K_{8}$, where local minima with low probabilities are increasingly
interspersed with other local minima corresponding to higher probabilities.

Figure 5: Scatter plots of probabilities of the correct Max-Cut solutions
against $\langle\hat{H}_{C}\rangle$ for individual minima of graphs $K_{6}$,
$K_{7}$ and $K_{8}$ for $L=1$ (in black), $L=2$ (in red) and $L=3$ (in green).

To further evaluate the performance of QAOA ansätze for various $L$ based on
the databases of minima and their respective Max-Cut probabilities and
$\langle\hat{H}_{C}\rangle$ values, we introduce the weighted average metric
$F$:

$F=\frac{1}{M|{\langle\hat{H}_{C}\rangle}_{min}|}{\sum_{m=1}^{M}}|{\langle\hat{H}_{C}\rangle}_{min}-{\langle\hat{H}_{C}\rangle}_{m}|[1-p_{m}(\ket{s})],$
(9)

where ${\langle\hat{H}_{C}\rangle}_{min}$ is the ${\langle\hat{H}_{C}\rangle}$
value of the global minimum. This formulation of $F$ is advantageous for two
reasons. First, it distinguishes circuit ansätze with $L_{ad}$ and $L_{min}$
layers, since for $L_{ad}$, $F=0$ because only the global minimum is present,
compared to $L_{min}$ where other local minima with lower Max-Cut
probabilities are present. Second, $F$ reflects the contribution of minima for
circuit ansätze with $L<L_{min}$ more accurately, since it is possible for
minima with lower ${\langle\hat{H}_{C}\rangle}$ values, including the global
minimum, to possess lower Max-Cut probabilities than their counterparts with
higher ${\langle\hat{H}_{C}\rangle}$ values. As $L$ increases, a decrease in
the value of $F$ can generally be interpreted as an improvement in QAOA
performance, since it corresponds to an increase in the proportion of local
minima with better probabilities. However, the converse situation, where $F$
increases as $L$ increases, may not necessarily signify a drop in QAOA
performance, as the well-funnelled organisation of the cost landscape, as well
as the guarantee of obtaining a better Max-Cut probability, may outweigh the
trade-off in obtaining a lower proportion of local minima with relatively good
probabilities. Nevertheless, we propose choosing circuit ansätze with $L$
layers that feature sufficiently low values of $F$ when simulating QAOA on
noisy quantum devices, as choosing circuits with higher $L$ may also increase
the impact of quantum gate and qubit decoherence noise.

Table 2:  The weighted average metric $F$ for graphs $K_{5}$ to $K_{8}$ of varying $L$. Graph | $L=1$ | $L=2$ | $L=3$
---|---|---|---
$K_{5}$ | 0.002276 | 0.067219 | 0.000000
$K_{6}$ | 0.083438 | 0.060651 | 0.007907
$K_{7}$ | 0.005269 | 0.119087 | 0.026351
$K_{8}$ | 0.081604 | 0.190701 | 0.038062

Analysing the $F$ values for $K_{5}$ to $K_{8}$ (Table II), we observe that
for $K_{5}$, $F=0$ for $L_{ad}=3$, differentiating it from $L_{min}=2$, as
expected. Interestingly, with the exception of $K_{6}$, $F$ appears to
increase for $L=2$ before decreasing for $L=3$. The increase in $F$ can mainly
be attributed to the general increase in the number of local minima for $L=2$
with comparably lower probabilities than that of the single global minimum for
$L=1$. This trend is evident in the scatter plots of the probabilities of the
correct Max-Cut solution against $\langle\hat{H}_{C}\rangle$ for the databases
of minima for varying $L$ (Fig. 5). The somewhat triangular convex hull of the
solution space, $C_{s}$, seems to become better defined with the transition
from $L=2$ to $L=3$. We also see a proliferation of local minima towards the
apex of the global minimum, which would explain the observed subsequent
decrease in $F$. Thus, a choice of $L=3$ would be adequate for solving the
Max-Cut problem for graphs $K_{5}$ to $K_{8}$ based on their $F$ values.

Figure 6:  3-regular unweighted graphs 6a–8e investigated in this study. Table 3:  Number of minima $M$ (top value), HCMPs (middle value) and $F$ (bottom value) for graphs 6a–8e of varying $L$ obtained from basin-hopping global optimisation. HCMPs with asterisks indicate that they correspond to the next highest local minimum rather than the global minimum. Graph | $L=1$ | $L=2$ | $L=3$ | $L=4$
---|---|---|---|---
6a | 1 | 5 | 23 | 145
0.400816 | 0.720917 | 0.933445* | 0.996304
0.254240 | 0.224891 | 0.167129 | 0.093993
6b | 3 | 20 | 191 | 1451
0.274835 | 0.659823 | 0.915348 | 0.978625
0.772861 | 0.615867 | 0.510817 | 0.340108
8a | 1 | 4 | 16 | 83
0.142701 | 0.349371 | 0.616672 | 0.748746
0.427166 | 0.363476 | 0.497637 | 0.428079
8b | 1 | 4 | 15 | 97
0.232056 | 0.420045* | 0.638057 | 0.767138
0.355057 | 0.325065 | 0.426828 | 0.380643
8c | 1 | 4 | 18 | 82
0.077352 | 0.182753 | 0.536762 | 0.701246
0.516636 | 0.560946 | 0.519856 | 0.429701
8d | 1 | 6 | 31 | 151
0.186302 | 0.520680 | 0.871573 | 0.972013
0.500506 | 0.360078 | 0.258910 | 0.141219
8e | 1 | 4 | 15 | 110
0.321737 | 0.574918* | 0.769975 | 0.918878*
0.286668 | 0.241691 | 0.302882 | 0.300042
Figure 7: Scatter plots of probabilities of the correct Max-Cut solutions
against $\langle\hat{H}_{C}\rangle$ for connected minima of graphs 6a–8e for
$L=1$ (in black), $L=2$ (in red), $L=3$ (in green) and $L=4$ (in blue).

###  III.2 3-regular unweighted graphs

Next, we analysed all connected 3-regular unweighted graphs with six and eight
vertices [66] from $L=1$ to $L=4$, labelled 6a–8e, respectively (Fig. 6). In
terms of the number of minima obtained from basin-hopping runs, the 3-regular
graphs generally possess lower $M$ values than the complete graphs of $K_{6}$
and $K_{8}$, particularly for graphs with eight vertices (Table III). For
graph 6b, starting from $L=1$ its QAOA ansatz gives rise to more than one
minimum, which subsequently produces a more rapid increase in $M$ compared to
its counterpart 6a. This phenomenon can largely be attributed to the
relatively more complex analytic expression of $\langle\hat{H}_{C}\rangle$
with the $L=1$ circuit ansatz [67] for 6b. Another interesting pair of graphs
is 8d and 8e, where, although their $L=1$ analytic expressions of
$\langle\hat{H}_{C}\rangle$ are identical and thus give rise to the same
$\langle\hat{H}_{C}\rangle_{min}$ value [67], we find that the trends of $M$
and HCMPs for higher $L$ are markedly different. Although the 3-regular series
of graphs give rise to lower $M$ for similar $L$, and hence somewhat simpler
energy landscapes compared to their complete graph counterparts, their HCMPs
are also comparatively lower, hence requiring a greater $L$ to achieve a
sufficiently high HCMP. In some cases, the HCMPs were derived not from the
global minimum, but from the next highest local minimum. (See Appendix B
detailing the differences in Max-Cut probabilities and
$\langle\hat{H}_{C}\rangle$ values between the two minima for these cases.)
This phenomenon appears sporadically without much predictability, most notably
for graph 8e, where the HCMP corresponds to the next highest local minimum for
$L=2$ and $L=4$. Overall, these observations further underline the importance
of evaluating QAOA performance based on the correct Max-Cut probabilities of
individual minima independently alongside their $\langle\hat{H}_{C}\rangle$
values.

The distributions of minima for the 3-regular graphs with varying $L$ differ
significantly from their complete graph counterparts (Fig. 7). We observe that
the convex hulls of the 3-regular graphs tend to take on a more compact shape,
with greater correlation between the $p(\ket{s})$ and
$\langle\hat{H}_{C}\rangle$ values for the individual minima. However, for
graphs 8a, 8b and 8e there is a notable absence of minima with intermediate
$p(\ket{s})$ and $\langle\hat{H}_{C}\rangle$ values, particularly for higher
$L$. This structure is also reflected in their disconnectivity graphs (refer
to Appendix C for the disconnectivity graphs of 6a-8e). Another major
difference of the 3-regular graphs is the much reduced energy differences
between minima with low proximity to the global minimum and their connected
transition states, producing more streamlined and single-funnelled
disconnectivity graphs than for $K_{6}$ and $K_{8}$. However, even though the
energy landscapes of the 3-regular graphs appear less complex and easier to
navigate than their complete graph counterparts, the local minima in their
energy landscapes give rise to a larger range of $p(\ket{s})$. Hence a greater
proportion of local minima with high energies possess low Max-Cut
probabilities. This trend is captured by comparing $F$ values between the
3-regular graphs and the complete graphs $K_{6}$ and $K_{8}$, where the former
graphs typically have much higher $F$ values than the latter.

Comparing $F$ values among the 3-regular graphs, we observe that graphs 6a and
6b follow a smooth downward trend with increasing $L$, while the 8-vertex
graphs tend to peak at $L=2$ and $L=3$ before decreasing, with the exception
of 8d, which follows a similar trend of the 6-vertex graphs. The performance
of the 8-vertex graphs at $L=2$ and $L=3$ can be attributed to an increase in
the number of local minima with low Max-Cut probabilities that outweigh the
general improvement in HCMPs and low-lying minima, while for graph 8d this
effect is reversed, with an increase in HCMPs and minima with good
probabilities. At $L=4$, a greater proportion of minima with high $p(\ket{s})$
appear for all 3-regular graphs, and we therefore recommend a minimum of $L=4$
when employing QAOA for these cases.

Table 4:  Expectation thresholds $d_{1}$ (top value) and $d_{2}$ (bottom value) for graphs 6a–8e of varying $L$. Graph | $L=2$ | $L=3$ | $L=4$
---|---|---|---
6a | 0.313189 | 0.546894 | 0.578616
0.643070 | 1.254700 | 1.342016
6b | 0.640787 | 1.453239 | 1.660319
0.951659 | 2.160089 | 2.470976
8a | - | 0.215689 | 0.338515
- | 0.673507 | 1.251349
8b | - | 0.250008 | 0.408135
- | 0.886500 | 1.652028
8c | - | 0.046924 | 0.315151
- | 0.183791 | 1.203264
8d | 0.052613 | 1.580115 | 1.857498
0.202901 | 2.457872 | 2.962362
8e | 0.096120 | 0.779749 | 0.645478
0.731083 | 1.578323 | 1.977320

Another factor that supports the choice of $L=4$ comes from the construction
of heuristic expectation thresholds that aim to identify minima with
sufficiently high $p(\ket{s})$ values. This analysis can be carried out by
finding the intercepts of the corresponding convex hulls with a suitable
probability cutoff $p_{op}$. For the 3-regular graphs we choose $p_{op}=0.5$
and define the difference in $\langle\hat{H}_{C}\rangle$ values from the
global minimum to the two intercepts as the worst-case and best-case
expectation cutoffs, $d_{1}$ and $d_{2}$, respectively, where $d_{1}<d_{2}$.
We observe that the expectation thresholds generally expand as $L$ increases
(Table IV), with $d_{1}$ and $d_{2}$ attaining their highest values at $L=4$.
For $L\leq{L_{min}}$ the widening and stabilising of expectation thresholds is
significant, along with the increase in $M$ as $L$ increases. We see that a
greater number of minima that possess a wider range of
$\langle\hat{H}_{C}\rangle$ values with a sufficiently high $p(\ket{s})$ exist
within the solution landscape for the QAOA ansatz.

$v_{1}$$v_{2}$$v_{4}$$v_{3}$3214$x$ Figure 8: Four-vertex weighted graph with a variable central weight $x$. The graphs $(G_{2},G_{3},G_{4},G_{5})$ correspond to $x=(0,3,4,5)$ respectively. Table 5:  Number of minima $M$ (top value), HCMPs (middle value) and $F$ (bottom value) for graphs $G_{2}-G_{5}$ with varying $L$ obtained from basin-hopping global optimisation. HCMPs with asterisks indicate that they were obtained from the next highest local minimum instead of the global minimum. Graph | $L=1$ | $L=2$ | $L=3$
---|---|---|---
$G_{2}$ | 3 | 16 | 87
0.477824 | 0.866404 | 0.981615
0.594432 | 0.338167 | 0.163886
$G_{3}$ | 7 | 109 | 1835
0.370522* | 0.784475 | 0.948704*
0.526539 | 0.300520 | 0.171897
$G_{4}$ | 4 | 36 | 308
0.540426* | 0.894074 | 1.000000
0.461355 | 0.245864 | 0.111340
$G_{5}$ | 9 | 183 | 3426
0.468114* | 0.881732* | 0.972620*
0.498733 | 0.328245 | 0.174754

###  III.3 Competing QAOA Max-Cut solutions

Finally, we explore competing Max-Cut solutions $\ket{\alpha\beta\alpha\beta}$
and $\ket{\alpha\alpha\beta\beta}$ for a series of four-vertex weighted graphs
with a common variable weight $x$, where $x=(0,3,4,5)$ correspond to the
graphs $(G_{2},G_{3},G_{4},G_{5})$, respectively (Fig. 8). These graphs open
up two modes of analysis: they allow comparison between $G_{2}$ with the more
complex graphs $G_{3}-G_{5}$, particularly with $G_{3}$, since both sets of
graphs have $\ket{s}=\ket{\alpha\beta\alpha\beta}$. Comparisons between
$G_{3}-G_{5}$ can also be carried out, since $G_{5}$ possesses a different
correct Max-Cut solution of $\ket{\alpha\alpha\beta\beta}$, while for $G_{4}$
both $\ket{\alpha\beta\alpha\beta}$ and $\ket{\alpha\alpha\beta\beta}$ are
correct Max-Cut solutions. We will denote the alternative Max-Cut solution for
a given graph as $\ket{t}$, thus for $G_{3}$,
$\ket{t}=\ket{\alpha\alpha\beta\beta}$; for $G_{5}$,
$\ket{t}=\ket{\alpha\beta\alpha\beta}$, and $G_{4}$ has no alternative
solution.

We find that implementing QAOA for the weighted graphs $G_{2}-G_{5}$ is more
difficult than for the complete unweighted graph $K_{4}$, as their energy
landscapes are much more complex due to an increase in $M$ and a decrease in
their respective HCMPs (Table V). The disconnectivity graphs of $G_{2}-G_{5}$
exhibit similar topological features to the 3-regular unweighted graphs,
possessing a well-funnelled organisation and minima featuring a wide range of
$\langle\hat{H}_{C}\rangle$ and Max-Cut probabilities (see Appendix D for the
disconnectivity graphs of $G_{2}-G_{5}$ coloured based on $p(\ket{s})$, and
Appendix E for the disconnectivity graphs of $G_{3}$ and $G_{5}$ coloured
based on $p(\ket{t})$). Unsurprisingly, $G_{2}$ has a lower $M$ and thus a
simpler landscape than the more strained graphs $G_{3}-G_{5}$, although its
collection of minima with modest Max-Cut probabilities produces a
comparatively high $F$ value up to $L=3$. In the range $G_{3}-G_{5}$, it is
interesting that even though $G_{5}$ with a different Max-Cut solution has a
more complex energy landscape than $G_{3}$, it yields a comparatively higher
HCMP, while $G_{4}$ exhibits the best QAOA performance as it has two distinct
Max-Cut solutions. The phenomenon where the HCMP arises for the next highest
local minimum rather than the global minimum was also observed for
$G_{3}-G_{5}$, especially for $G_{5}$ (refer to Appendix B for differences in
Max-Cut probabilities and $\langle\hat{H}_{C}\rangle$ values between the two
minima for these cases). All four graphs feature a monotonic decrease in $F$
as $L$ increases. Hence, as for $K_{5}-K_{8}$, we recommend choosing $L=3$ in
solving the Max-Cut problem for $G_{2}-G_{5}$. Finally, it is noteworthy that
for $G_{4}$ and the $L=4$ ansatz, the number of minima greatly exceeds that of
the $L_{min}=3$ ansatz, in contrast to the behaviour observed for $K_{5}$,
highlighting the general unpredictability of the energy landscape complexity
after $L_{min}$.

Table 6:  Expectation thresholds $d_{1}$ (top value) and $d_{2}$ (bottom value) for graphs $G_{3}-G_{5}$ of varying $L$. Graph | $L=2$ | $L=3$
---|---|---
$G_{3}$ | 0.099242 | 0.433945
1.124694 | 2.008400
$G_{4}$ | 0.489057 | 0.536238
1.447505 | 1.722258
$G_{5}$ | 0.217591 | 0.459177
1.566966 | 1.896063
Figure 9: Scatter plots of probabilities of competing Max-Cut solutions
$\ket{\alpha\beta\alpha\beta}$ (in red) and $\ket{\alpha\alpha\beta\beta}$ (in
blue) against $\langle\hat{H}_{C}\rangle$ for connected minima of graphs
$G_{3}$, $G_{4}$ and $G_{5}$ of varying $L$.

For graphs $G_{3}-G_{5}$, the scatter plots of the $\langle\hat{H}_{C}\rangle$
values for connected minima and Max-Cut probabilities were used to construct
both the convex hulls of the correct Max-Cut solution $C_{s}$ and the
alternative solution $C_{t}$ (Fig. 9). We find that, similar to the unweighted
graphs, both $C_{s}$ and $C_{t}$ take on more definite shapes as they become
populated with more minima for increasing $L$. This trend allowed us to
investigate $d_{1}$ and $d_{2}$ by finding the intercepts of $C_{s}$ with the
horizontal line $p_{op}$, setting $p_{op}=0.5$ for $G_{3}$ and $G_{5}$, and
$p_{op}=0.25$ for $G_{4}$ (Table VI). Additionally, one may also use the
intercepts between $C_{s}$ and the left edge of $C_{t}$ as expectation and
threshold cut-offs to identify higher-quality minima with a high probability
of finding $\ket{s}$ and a low probability of obtaining $\ket{t}$. (See
Appendix F for more details.) As for the 3-regular graphs, the divergence of
$d_{1}$ and $d_{2}$ with increasing $L$ makes $L=3$ a good choice for graphs
$G_{3}-G_{5}$.

##  IV Conclusion

In this work we explore the solution landscapes of QAOA ansätze applied to a
variety of weighted and unweighted graphs by means of the energy landscapes
framework, using disconnectivity graphs to visualise their topological
features. We find that the corresponding landscapes are largely funnelled,
suggesting that location of low-lying minima should not be particularly
difficult. Under practical conditions when simulating QAOA on a quantum
device, the optimisation regime is thus more likely to find a minimum close to
the global minimum with a good correct Max-Cut solution probability. Even
under the worst-case scenario where each experiment finds a different local
minimum, so long as the local minimum is sufficiently close to the global
minimum, a significant proportion of the number of shots per experiment will
correspond to the Max-Cut solution. This result further demonstrates the
robustness of QAOA in solving the Max-Cut problem.

We have also developed a weighted average metric $F$ to evaluate the
performance of QAOA ansätze from their corresponding databases of minima. This
parameter allows one to choose a suitable number of circuit layers that
balances the likelihood of obtaining good solution probabilities from local
minima with an adequate circuit depth that minimises the impact from quantum
noise.

Finally, we have established two ways in which expectation thresholds can be
established to determine the cut-off for minima with high $p(\ket{s})$. The
solution landscapes we have characterised suggest that QAOA is a good VQA
candidate to demonstrate practical quantum advantage. In future work we plan
to extend these results to quantum machine learning (QML) algorithms, such as
variational quantum classifiers (VQCs), which minimise a given cost function
to classify data [68].

## Software Availability

The GMIN, OPTIM and PATHSAMPLE programs are available for use under the Gnu
General Public License. They can be downloaded from the Cambridge Landscape
Database at www-wales.ch.cam.ac.uk.

## Appendix A: Basin-hopping global optimisation with GMIN

For each basin-hopping optimisation run, we performed 10,000 basin-hopping
steps for the unweighted graphs $K_{3}-K_{8}$, 6a$-$8e, and for the weighted
graphs $G_{2}-G_{5}$, using the GMIN program for varying $L$. Each local
minimisation had a minimum root-mean squared (RMS) gradient convergence
criterion of $1.0\times 10^{-10}$ a.u., where the analytic gradients of the
parameterised rotation gates of the cost and mixer layers of the QAOA ansatz
were evaluated with the parameter-shift rule using Eq. 7. To accept/reject
basin-hopping steps we employed a Metropolis criterion with a basin-hopping
temperature of 1.0 a.u. If the minimum at step $j$ has an expectation value
$\langle\hat{H}_{C}\rangle_{j}$ that is lower than the preceding iteration,
i.e. $\langle\hat{H}_{C}\rangle_{j}<\langle\hat{H}_{C}\rangle_{j-1}$,
$\langle\hat{H}_{C}\rangle_{j}$ then the corresponding angular coordinates
$\bm{\theta}_{j}$ are accepted and used for the next step. If
$\langle\hat{H}_{C}\rangle_{j}\geq\langle\hat{H}_{C}\rangle_{j-1}$, then
$\langle\hat{H}_{C}\rangle_{j}$ and $\bm{\theta}_{j}$ is accepted with a
probability of
$\exp({-(\langle\hat{H}_{C}\rangle_{j}-\langle\hat{H}_{C}\rangle_{j-1})/{k}T})$.
Otherwise, the new minimum is rejected.

Basin-hopping moves were proposed by random perturbations of up to 1.0 rad for
each angular coordinate in $\bm{\theta}_{j}$. At the end of each run, the
collection of minima that differ by at least $1.0\times 10^{-9}$ a.u in their
$\langle\hat{H}_{C}\rangle$ values were saved to provide a starting database
for construction of the energy landscape using the OPTIM and PATHSAMPLE
programs.

## Appendix B: Differences in Max-Cut probability and expectation values for
non-global HCMP cases

Table 7:  Expectation differences $\Delta\langle\hat{H}_{C}\rangle$ and Max-Cut probability differences ${\Delta}p(\ket{s})$ between the global and next highest local minima for non-global HCMP cases. Graph | $L$ | $\Delta\langle\hat{H}_{C}\rangle$ | ${\Delta}p(\ket{s})$
---|---|---|---
6a | 3 | 0.037856 | 0.010121
8b | 2 | 0.075696 | 0.040029
8e | 2 | 0.046450 | 0.052876
8e | 4 | 0.006337 | 0.009846
$G_{3}$ | 1 | 0.127864 | 0.096306
$G_{3}$ | 3 | 0.015904 | 0.003739
$G_{4}$ | 1 | 1.142149 | 0.076922
$G_{5}$ | 1 | 0.342188 | 0.231787
$G_{5}$ | 2 | 0.036357 | 0.135824
$G_{5}$ | 3 | 0.001111 | 0.001329

Table VII shows the expectation differences
$\Delta\langle\hat{H}_{C}\rangle=\langle\hat{H}_{C}\rangle_{min}-\langle\hat{H}_{C}\rangle_{\text{HCMP}}$
and correct Max-Cut probability differences $\Delta
p(\ket{s})=\text{HCMP}-p_{\text{GM}}(\ket{s})$ for cases where the next
highest local minimum has a Max-Cut probability greater than the global
minimum. In general, cases with a large circuit depth $L$ tend to have lower
$\Delta\langle\hat{H}_{C}\rangle$ and $\Delta p(\ket{s})$ as their HCMPs
approach the optimal value of 1. As mentioned in the main text this phenomenon
appears to occur sporadically, as exemplified by 8e and $G_{3}$, where it does
not occur for the intermediate circuit depth of $L=3$. A notable case is
$G_{5}$, which features a higher Max-Cut probability in the next highest local
minimum for all sampled $L$. We propose to investigate this phenomenon in
future work to see if it arises systematically for particular classes of
connected graphs.

## Appendix C: Disconnectivity graphs of the 3-regular graphs

Figure 10: Disconnectivity graphs for 6a with varying circuit depth $L$.

Figure 11: Disconnectivity graphs for 6b with varying circuit depth $L$.

Figure 12: Disconnectivity graphs for 8a with varying circuit depth $L$.

Figure 13: Disconnectivity graphs for 8b with varying circuit depth $L$.

Figure 14: Disconnectivity graphs for 8c with varying circuit depth $L$.

Figure 15: Disconnectivity graphs for 8d with varying circuit depth $L$.
Figure 16: Disconnectivity graphs for 8e with varying circuit depth $L$.

## Appendix D: Disconnectivity graphs of the correct Max-Cut solutions for
$G_{2}-G_{5}$

Figure 17: Disconnectivity graphs for $G_{2}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the optimal Max-Cut
state of $\ket{\alpha\beta\alpha\beta}$.

Figure 18: Disconnectivity graphs for $G_{3}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the optimal Max-Cut
state of $\ket{\alpha\beta\alpha\beta}$.

Figure 19: Disconnectivity graphs for $G_{4}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the optimal Max-Cut
states of $\ket{\alpha\beta\alpha\beta}$ and $\ket{\alpha\alpha\beta\beta}$.

Figure 20: Disconnectivity graphs for $G_{5}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the optimal Max-Cut
state of $\ket{\alpha\alpha\beta\beta}$.

## Appendix E: Disconnectivity graphs of the alternative Max-Cut solutions for
$G_{3}$ and $G_{5}$

Figure 21: Disconnectivity graphs for $G_{3}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the opposing Max-Cut
state of $\ket{\alpha\alpha\beta\beta}$.

Figure 22: Disconnectivity graphs for $G_{5}$ with varying circuit depth $L$.
Minima are coloured based on the probability of obtaining the opposing Max-Cut
state of $\ket{\alpha\beta\alpha\beta}$.

## Appendix F: Expectation and probability thresholds from the convex hull
intercepts for $G_{3}$ and $G_{5}$

Table 8:  Expectation thresholds $d_{3}$ (top value) and $d_{4}$ (bottom value) for graphs $G_{3}$ and $G_{5}$ of varying $L$, where $d_{3}<d_{4}$. Graph | $L=2$ | $L=3$
---|---|---
$G_{3}$ | 0.111125 | 0.439036
0.546385 | 0.794808
$G_{5}$ | 0.329259 | 0.476747
0.856478 | 0.819369
Table 9:  Max-Cut probability thresholds $p_{1}$ (top value) and $p_{2}$ (bottom value) for graphs $G_{3}$ and $G_{5}$ of varying $L$, corresponding to the expectation thresholds $d_{3}$ and $d_{4}$ respectively. Graph | $L=2$ | $L=3$
---|---|---
$G_{3}$ | 0.465937 | 0.495366
0.665681 | 0.799744
$G_{5}$ | 0.373801 | 0.482335
0.680284 | 0.790380

Tables VIII and IX summarise the two intercepts
$(\langle\hat{H}_{C}\rangle_{min}+d_{3},p_{1})$ and
$(\langle\hat{H}_{C}\rangle_{min}+d_{4},p_{2})$ of $C_{s}$ with the left edge
of the triangular convex hull $C_{t}$ for graphs $G_{3}$ and $G_{5}$, as shown
in Fig. 9. This summary provides an alternative method to establish
expectation cut-offs if additional information for alternative Max-Cut states
is available. As for the expectation thresholds $d_{1}$ and $d_{2}$ in Table
VI, $d_{3}$ and $d_{4}$ also increase with increasing $L$, which further
justifies the choice of $L=3$ in refining minima within $C_{s}$ that possess
both good correct Max-Cut probabilities and low opposing Max-Cut
probabilities. The probability cut-off of $p_{op}=0.5$ used in Table VI also
appears to be situated between the probability thresholds $p_{1}$ and $p_{2}$,
suggesting that this choice of $p_{op}$ is optimal.

## References

  * O’Malley  _et al._ [2016] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, and N. D. et al., Scalable quantum simulation of molecular energies, Phys. Rev. X. 6, 031007 (2016).
  * Leymann and Barzen [2020] F. Leymann and J. Barzen, The bitter truth about gate–based quantum algorithms in the NISQ era, Quant. Sci. Tech. 5, 044007 (2020).
  * Skosana and Tame [2021] U. Skosana and M. Tame, Demonstration of Shor’s factoring algorithm for ${N}=21$ on IBM quantum processors, Sci. Rep. 11, 16599 (2021).
  * Burnett  _et al._ [2019] J. J. Burnett, A. Bengtsson, M. Scigliuzzo, D. Niepce, M. Kudra, P. Delsing, and J. Bylander, Decoherence benchmarking of superconducting qubits, npj Quantum Inf. 5, 16599 (2019).
  * Chen  _et al._ [2023] S. Chen, J. Cotler, H. Huang, and J. Li, The complexity of NISQ, Nat. Commun. 14, 6001 (2023).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru–Guzik, The theory of variational hybrid quantum–classical algorithms, New J. Phys. 18, 023023 (2016).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M. Yung, X. Zhou, P. J. Love, A. Aspuru–Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 4213 (2014).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera–Lierta, T. H. Kyaw, T. Haug, S. Alperin–Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, and T. M. et al., Noisy intermediate–scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019).
  * Ravi  _et al._ [2022] G. Ravi, K. N. Smith, P. Gokhale, A. Mari, N. Earnest, A. Javadi–Abhari, and F. T. Chong, VAQEM: A variational approach to quantum error mitigation, IEEE , 288–303 (2022).
  * Ito  _et al._ [2023] K. Ito, W. Mizukami, and K. Fujii, Universal noise–precision relations in variational quantum algorithms, Phys. Rev. Res. 5, 023025 (2023).
  * Fontana  _et al._ [2021] E. Fontana, N. Fitzpatrick, D. M. Ramo, R. Duncan, and I. Rungger, Evaluating the noise resilience of variational quantum algorithms, Phys. Rev. A 104, 022403 (2021).
  * Halder  _et al._ [2023] S. Halder, C. Shrikhande, and R. Maitra, Development of zero–noise extrapolated projective quantum algorithm for accurate evaluation of molecular energetics in noisy quantum devices, J. Chem. Phys. 159, 114115 (2023).
  * Weidinger  _et al._ [2023] A. Weidinger, G. B. Mbeng, and W. Lechner, Error mitigation for quantum approximate optimization, Phys. Rev. A 108, 032408 (2023).
  * Bravyi  _et al._ [2022] S. Bravyi, O. Dial, J. M. Gambetta, D. Gil, and Z. Nazario, The future of quantum computing with superconducting qubits, J. Appl. Phys. 132, 160902 (2022).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, and K. T. et al., Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500–505 (2023).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014).
  * Zhang  _et al._ [2022] Y. Zhang, L. Cincio, C. F. A. Negre, P. Czarnik, P. J. Colesa, P. M. Anisimov, S. M. Mniszewski, S. Tretiak, and P. A. Dub, Variational quantum eigensolver with reduced circuit complexity, npj Quantum Inf. 8 (2022).
  * Wurtz and Love [2022] J. Wurtz and P. J. Love, Counterdiabaticity and the quantum approximate optimization algorithm, Quantum 6, 635 (2022).
  * Buonaiuto  _et al._ [2023] G. Buonaiuto, F. Gargiulo, G. D. Pietro, M. Esposito, and M. Pota, Best practices for portfolio optimization by quantum computing, experimented on real quantum devices, Sci. Rep. 13 (2023).
  * Kremenetski  _et al._ [2021] V. Kremenetski, T. Hogg, S. Hadfield, S. J. Cotton, and N. M. Tubman, Quantum alternating operator ansatz (QAOA) phase diagrams and applications for quantum chemistry (2021).
  * Hadfield  _et al._ [2022] S. Hadfield, T. Hogg, and E. G. Rieffel, Analytical framework for quantum alternating operator ansätze, Quantum Sci. Tech. 8, 015017 (2022).
  * Wang  _et al._ [2021] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise–induced barren plateaus in variational quantum algorithms, Nat. Commun. 12, 6961 (2021).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nat. Commun. 9, 4812 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nat. Commun. 12, 1791 (2021).
  * Wierichs  _et al._ [2020] D. Wierichs, C. Gogolin, and M. Kastoryano, Avoiding local minima in variational quantum eigensolvers with the natural gradient optimizer, Phys. Rev. Res. 2, 043246 (2020).
  * Anschuetz and Kiani [2022] E. R. Anschuetz and B. T. Kiani, Quantum variational algorithms are swamped with traps, Nat. Commun. 13 (2022).
  * Sack  _et al._ [2023] S. H. Sack, R. A. Medina, R. Kueng, and M. Serbyn, Recursive greedy initialization of the quantum approximate optimization algorithm with guaranteed improvement, Phys. Rev. A 107, 062404 (2023).
  * Anschuetz and Kiani [2018] E. R. Anschuetz and B. T. Kiani, Performance of the quantum approximate optimization algorithm on the maximum cut problem (2018).
  * Wales [2003] D. J. Wales,  _Energy landscapes_ (Cambridge University Press, Cambridge, 2003).
  * Wales [2018] D. J. Wales, Exploring energy landscapes, Ann. Rev. Phys. Chem. 69, 401 (2018).
  * Li and Scheraga [1987] Z. Li and H. A. Scheraga, Monte Carlo-minimization approach to the multiple-minima problem in protein folding, Proc. Natl. Acad. Sci. USA 84, 6611 (1987).
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms, J. Phys. Chem. A 101, 5111 (1997).
  * Wales and Scheraga [1999] D. J. Wales and H. A. Scheraga, Global optimization of clusters, crystals and biomolecules, Science 285, 1368 (1999).
  * Wales [2002] D. J. Wales, Discrete path sampling, Mol. Phys. 100, 3285 (2002).
  * Wales [2004] D. J. Wales, Some further applications of discrete path sampling to cluster isomerization, Mol. Phys. 102, 891 (2004).
  * Choy and Wales [2023] B. Choy and D. J. Wales, Molecular energy landscapes of hardware–efficient ansätze in quantum computing, J. Chem. Theory. Comput. 19, 1197–1206 (2023).
  * Burton  _et al._ [2023] H. G. A. Burton, D. Marti–Dafcik, D. P. Tew, and D. J. Wales, Exact electronic states with shallow quantum circuits from global optimisation, npj Quantum Inf. 9 (2023).
  * Wales  _et al._ [1998] D. J. Wales, M. A. Miller, and T. R. Walsh, Archetypal energy landscapes, Nature 394, 758–760 (1998).
  * Bogdan and Wales [2004] T. V. Bogdan and D. J. Wales, New results for phase transitions from catastrophe theory, J. Chem. Phys. 120, 11090 (2004).
  * D. J. Wales [2023] D. J. Wales, GMIN: A program for finding global minima and calculating thermodynamic properties from basin–sampling, https://www-wales.ch.cam.ac.uk/GMIN/ (2023), accessed 1 December 2023.
  * Mari  _et al._ [2021] A. Mari, T. R. Bromley, and N. Killoran, Estimating the gradient and higher–order derivatives on quantum hardware, Phys. Rev. A 103, 012405 (2021).
  * Broyden [1970] C. G. Broyden, The convergence of a class of double-rank minimization algorithms 1. General considerations, J. Inst. Math. Appl. 6, 76 (1970).
  * Fletcher [1970] R. Fletcher, A new approach to variable metric algorithms, Comput. J. 13, 317 (1970).
  * Goldfarb [1970] D. Goldfarb, A family of variable metric updates derived by variational means, Math. Comput. 24, 23 (1970).
  * Shanno [1970] D. F. Shanno, Conditioning of quasi-Newton methods for function minimization, Math. Comput. 24, 647 (1970).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory bfgs method for large scale optimization, Math. Program. 45, 503 (1989).
  * Nocedal [1980] J. Nocedal, Updating quasi-Newton matrices with limited storage, Math. Comput. 35, 773 (1980).
  * Paleico and Behler [2020] M. L. Paleico and J. Behler, A flexible and adaptive grid algorithm for global optimization utilizing basin hopping Monte Carlo, J. Chem. Phys. 152, 094109 (2020).
  * Noé and Fischer [2008] F. Noé and S. Fischer, Transition networks for modeling the kinetics of conformational change in macromolecules, Curr. Opin. Struc. Biol. 18, 154 (2008).
  * Prada-Gracia  _et al._ [2009] D. Prada-Gracia, J. Gómez-Gardenes, P. Echenique, and F. Falo, Exploring the free energy landscape: From dynamics to networks and back, PLoS Comput. Biol. 5, e1000415 (2009).
  * Wales [2010] D. J. Wales, Energy landscapes: some new horizons, Curr. Opin. Struct. Biol. 20, 3 (2010).
  * Trygubenko and Wales [2004] S. A. Trygubenko and D. J. Wales, A doubly nudged elastic band method for finding transition states, J. Chem. Phys. 120, 2082 (2004).
  * Sheppard  _et al._ [2008] D. Sheppard, R. Terrell, and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys. 128, 134106 (2008).
  * Mills  _et al._ [1995] G. Mills, H. Jónsson, and G. K. Schenter, Reversible work transition state theory: application to dissociative adsorption of hydrogen, Surf. Sci. 324, 305 (1995).
  * Jónsson  _et al._ [1998] H. Jónsson, G. Mills, and K. W. Jacobsen, Nudged elastic band method for finding minimum energy paths of transitions, in  _Classical and quantum dynamics in condensed phase simulations_ (World Scientific, 1998) Chap. 16, pp. 385–404.
  * Henkelman  _et al._ [2000] G. Henkelman, B. P. Uberuaga, and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113, 9901 (2000).
  * Henkelman and Jónsson [2000] G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113, 9978 (2000).
  * Munro and Wales [1999] L. J. Munro and D. J. Wales, Defect migration in crystalline silicon, Phys. Rev. B 59, 3969–3980 (1999).
  * Henkelman and Jónsson [1999] G. Henkelman and H. Jónsson, A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives, J. Chem. Phys. 111, 7010 (1999).
  * Kumeda  _et al._ [2001] Y. Kumeda, L. J. Munro, and D. J. Wales, Transition states and rearrangement mechanisms from hybrid eigenvector-following and density functional theory. application to C10H10 and defect migration in crystalline silicon, Chem. Phys. Lett. 341, 185 (2001).
  * Carr  _et al._ [2005] J. M. Carr, S. A. Trygubenko, and D. J. Wales, Finding pathways between distant local minima, J. Chem. Phys. 122, 234903 (2005).
  * Dijkstra [1959] E. W. Dijkstra, A note on two problems in connexion with graphs, Numerische Math. 1, 269 (1959).
  * Becker and Karplus [1997] O. M. Becker and M. Karplus, The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics, J. Chem. Phys. 106, 1495 (1997).
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, and S. B. et al., Quantum approximate optimization of non–planar graph problems on a planar superconducting processor, Nat. Phys. 17, 332–336 (2021).
  * Guerreschi and Matsuura [2019] G. G. Guerreschi and A. Y. Matsuura, QAOA for max–cut requires hundreds of qubits for quantum speed–up, Sci. Rep. 9, 332–336 (2019).
  * Fernández–Pendás  _et al._ [2022] M. Fernández–Pendás, E. F. Combarro, S. Vallecorsa, J. Ranilla, and I. F. Rúa, A study of the performance of classical minimizers in the Quantum Approximate Optimization Algorithm, J. Comput. Appl. Math. 404, 113388 (2022).
  * Jäger and Krems [2023] J. Jäger and R. V. Krems, Universal expressiveness of variational quantum classifiers and quantum kernels for support vector machines, Nat. Commun. 14, 576 (2023).


# Sr2VO3FeAs as Compared with Other Fe-based Superconductors

I.I. Mazin  Code 6393, Naval Research Laboratory, Washington, D.C. 20375

(Printed on )

###### Abstract

One of the most popular scenarios for the superconductivity in Fe-based
superconductors (FeBSC) posits that the bosons responsible for electronic
pairing are spin-fluctuations with a wave vector spanning the hole Fermi
surfaces (FSs) near $\Gamma$ and the electron FSs near M points. So far all
FeBSC for which neutron data are available do demonstrate such excitations,
and the band structure calculations so far were finding quasi-nested FSs in
all FeBSC, providing for a peak in the spin susceptibility at the desired wave
vectors. However, the newest addition to the family, Sr2VO3FeAs, has been
calculated to have a very complex FS with no visible quasi-nesting features.
It was argued therefore that this material does not fall under the existing
paradigm and calls for revisiting our current ideas about what is the likely
cause of superconductivity in FeBSC. In this paper, I show that the visible
complexity of the FS is entirely due to the V-derived electronic states.
Assuming that superconductivity in Sr2VO3FeAs, as in the other FeBSC,
originates in the FeAs layers, and the superconducting electrons are sensitive
to the susceptibility of the FeAs electronic subsystem, I recalculate the bare
susceptibility, weighting the electronic states with their Fe character, and
obtain a susceptibility that fully supports the existing quasi-nesting model.
Besides, I find that the mean-filed magnetic ground state is the checkerboard
in V sublattuce and stripes in the Fe sublattice.

The recently discovereddisc  Fe-based high-temperature superconductors (FeBSC)
represent a challenging case for the theory of superconductivity. They appear
to be rather different from cuprates in terms of their electronic structure,
magnetic order, correlation effects, and superconducting symmetryreview . So
far the most popular suggestion for the pairing mechanism has been one that
assigns the role of an intermediate boson to spin fluctuations with wave
vectors close to Q=($\pi,\pi)$ (in the two-Fe Brillouin zone). There are two
ways to generate such spin fluctuations: one assumes superexchange between the
second neighbors in the Fe lattice and the other exploits the fact that the
non-interacting spin susceptibility calculated using the one-electron band
structure has a peak, or better to say a broad maximum close to ($\pi,\pi)$
(see review Ref. review, ). A strong argument in favor of the latter scenario
was the case of FeSe, where the parent magnetic compound FeTe shows an
antiferromagnetic order at a different wave vector. both in the experiment and
in the calculations, but the calculated spin susceptibility is still peaked
Q=($\pi,\pi),$ and the experiment also observes spin fluctuations with the
same wave vector. Also, the fact that FeBSC lack strong Coulomb
correlationsAnis ; Tom  speaks against the former alternative.

Recently, however, a new FeBSC, Sr2VO3FeAs, has been discovered which
seemingly violates this so far meticulously observed rule. The calculated
Fermi surface (FS)WEP  appears to be much more complex than in the other
investigated FeBSC, and there is no visual indication of any quasinesting
topology. Lee and PickettWEP  argued that Sr2VO3FeAs represents “a new
paradigm for Fe-pnictide superconductors”, and inferred that “there is no
reason to expect an s± symmetry of superconducting order parameter ($i.e.$ a
different sign on the two FSs) in Sr2VO3FeAs.

Figure 1: The Fermi surfaces of Sr2VO3FeAs. The $\Gamma$ points are in the
corners, the M point in the center of the shown Brillouin zone. The colored
(dark) portion are the parts with the predominantly Fe character. The rest is
predominantly V. (color online)

I have repeated the calculations of Lee and Pickett and have obtained the FS
that was similar to theirsnote  (Fig. 1). I have also verified that the bare
susceptibility without any account for the matrix elements

$\chi_{0}(\mathbf{q)}=-\sum_{\mathbf{k\alpha\beta}}\frac{f(\varepsilon_{\mathbf{k\alpha}})-f(\varepsilon_{\mathbf{k+q,\beta}})}{\varepsilon_{\mathbf{k\alpha}}-\varepsilon_{\mathbf{k+q,\beta}}+i\delta}$
(1)

indeed does not have any peak at Q=($\pi,\pi)$ (Fig. 2). In fact, it has a
peak at an entirely different wave vector, $(\pi,0.4\pi),$ as anticipated by
Lee and Pickett. However, this does not take into account the fact that the
calculated Fermi surface is really a superposition of two FS systems, one
originating from the FeAs planes, and the other from VO ones. While there is
some hybridization between the two systems of bands (at least along the XM
direction; see Ref. WEP,  for details), as well as a magnetic coupling and a
magnetic moment on V, and maybe even Coulomb correlation effects on V site,
electrons derived from the Fe $d$-orbitals couple mostly with the spin
fluctuations on the Fe sites. This is a simple consequence of the Hund’s rule.
With that in mind, I colored the parts of the Fermi surface in Fig. 1 that
have predominantly Fe character.

Figure 2: The bare susceptibility (the real part) calculated with a constant
matrix element independently of the wave function character. The band
structure had been averaged over $k_{z}$ before the integration. The corners
of the plot correspond to $\mathbf{q}=(0,0)$, $(\pi,0)$, $(0,\pi)$, and
$(\pi,\pi)$. The vertical scale is in arbitrary units. (color online)

Figure 3: The bare susceptibility calculated as in Fig.2, but with matrix
elements taken as the product of the Fe weights for the corresponding wave
functions. The top panel shows the real part, the bottom one the imaginary
part. (color online)

Imagine now that the unpainted parts of the FS disappear. What remains after
this mental tour de force closely resembles the familiar FSs of other FeBSC.
Taking into account the above argument regarding the special role of the Fe
spin fluctuations, we can rewrite Eq. 1 as

$\tilde{\chi}_{0}(\mathbf{q)}=-\sum_{\mathbf{k\alpha\beta}}\frac{f(\varepsilon_{\mathbf{k\alpha}})-f(\varepsilon_{\mathbf{k+q,\beta}})}{\varepsilon_{\mathbf{k\alpha}}-\varepsilon_{\mathbf{k+q,\beta}}+i\delta}A_{\mathbf{k\alpha}}A_{\mathbf{k+q,\beta}},$
(2)

where $A_{\mathbf{k\alpha}}$ is the relative weight of the Fe orbitals in the
$|\mathbf{k}\alpha\mathbf{>}$ wave function. The result (Fig. 3), as expected,
shows the same structure as for the other pnictides, especially for the real
part of susceptibility, which is the one relevant for superconductivity.

I conclude that, unfortunately, Sr2VO3FeAs, despite being an interesting and
in many aspects unusual FeBSC, does not represent a new paradigm, but rather
falls into the same class as other pnictides. It is also worth noting that
while it has been established both experimentallyAnis ; Tom  and
computationallyAnis ; Antoin  that the FeAs subsystem is only weakly
correlated, this had not been obvious a priori, and it is not obvious for the
V-O subsystem in Sr2VO3FeAs. Being essentially in a vanadium oxide layer (and
vanadium oxide is strongly correlated in the bulk form), V in Sr2VO3FeAs may
be subject to strong Hubbard correlations that would remove V states from the
Fermi levelLDAU . Thus, strictly speaking, the conclusion above should be
formulated as follows: even if Sr2VO3FeAs is a weakly correlated metal and the
FS calculated within the density functional theory is realistic, the fact that
the overall topology seems on the first glance to be different from other
pnictides is misleading and the spin fluctuation spectrum is likely to be
rather similar.

At the end, let me briefly touched upon a separate, but equally (if not more)
interesting issue of the magnetic ground state and magnetic properties of
Sr2VO3FeAs within the density functional theory (DFT). It is well knowreview
that DFT seriously overestimates the tendency to magnetism in FeBSCs, so that
the calculated ground state appears strongly antiferromagnetic even in the
materials that sho no long range magnetic order (phosphates, selenide). This
is routinely ascribed to the mean-filead character of DFT. However, it is of
course interesting to see what is the (magnetic) ground state in the mean
filed, even when in real life the ground state is paramagnetic. For all FeBSCs
studied so far the antiferromagnetic stripe magnetic structure is by far the
lowest in energy (energy gain of the order of 200 meV per Fe compared to a
nonmagnetic solution), while the ferromagnetic structure is barely stable if
at all.

Most likely, the DFT ground state of FeBSCs is also antiferromagnetic in-
plane. However, even the nonmagetic unit cell contains 16 atoms, which makes
it extremely difficult to investigate the energy landscape for possible
antiferromagnetic pattern. Thus, it makes sense to study possible
ferro(ferri)magnetic solutions, in hope to extract at least some useful
information. This approach was adapted in Ref. Shein  (although these authors
do not present any nonmagnetic calculations, actually relevant for
superconductivity). They found a solution with a moment on V ($\sim 1.5$
$\mu_{B}),$ but not on Fe. Lee and Pickett found another, ferrimagnetic
solution, with opposite moments on V and Fe, the former being largerLP . Using
different starting configurations, I was able to converge to three different
ground states within the same symmetry, as shown in the Table, as well as to
two lower-symmetry states, as illustrated in Fig. 4(b,c,d): interlayer
antiferromagnetic V sublattice, where the V layers are ferromagnetic, and
antiferromagnetically stacked, while Fe is nonmagnetic, and Fe-checkerboard,
where Fe forms a Neel plane and V in nonmagnetic. After that, I have
calculated two configurations in the double (four formula units) cell, which I
feel are the most relevant because of the superexchange interaction in the V
layers: V-checkerboard with nonmagnetic Fe, and V-checkerboard combined with
the stripe order in the Fe layers (Fig. 4)

Figure 4: Magnetic configurations used in the Table 1. Hollow symbols indicate
nonmagnetic atoms, blue (dark) spin-up moments, red (light) spin-down moments.
Circles are Fe atoms, upward and downward pointing triangles are two V layers
in the unit cell. The configurations are: (a) NM:nonmagnetic, (b)
FM:ferromagnetic, (c) half-FM, (d) FiM: ferrimagnetic (Fe and V spins are
antiparallel), (e) V-AF: antiferromagnetically stacked FM V layers,
nonmagnetic Fe (f) Fe-cb: checkerboard Fe planes, weakly ferromagnetic V
planes, (g) V-cb: checkerboard V planes, ferromagnetic Fe planes (h) V-cb
combined with Fe stripes. Minimal crystallographic unit cell is shown in each
case, and in the last panel dashed lines connect V atoms in the same layer
(color online)

A few observations are in place: (1) the state found in Ref. Shein  is not the
ground state even within that symmetry; (2) unlike all other FeBSCs, FeAs
planes can support a very stable ferromagnetic state; (3) the interaction
between V and Fe is ferromagnetic, that is, not of superexchange character,
(4) the magnetic coupling between V and Fe is so weak that V does not induce
any magnetization on Fe, unless one already starts with a magnetic Fe; (5) It
is more important, from the total energy point of view, to have magnetic
moment on V that on Fe — a bit surprising, given that V has a weaker Hund’s
rule coupling; (6) V sublattice itself has a net antiferromagnetic
interaction: if Fe is not magnetic, V orders antiferromagnetically; (7) Unless
some more exotic ground state will be discovered, the total energy is
minimized when V layers order in the Neel (checkerboard) fashion, while Fe
orders the same way as in other pnictides, forming stripes; (8) most
importantly, a number of very different magnetic states are nearly degenerate
in energy. This last fact may be the key to the experimental fact that the
actual material is paramagnetic despite the fact that on the mean field level
it is more magnetic than other pnictides. This is an extremely intriguing
situation and the magnetism Sr2VO3FeAs deserves a more elaborated experimental
and theoretical study that is beyond the scope of this paper.

Table 1: Properties of some stable magnetic solutions in the Generalized Gradient Approximation of the DFT. All energies are given with respect to the nonmagnetic state. The magnetic states are described in Fig. 4. For the V-cb configuration I was able to converge to two different solution, high-spin and low-spin, with essentially the same energies |  $M_{Fe},$ $\mu_{B}$ |  $M_{V},$ $\mu_{B}$ |  $\Delta E.$ meV/Fe
---|---|---|---
FM | 2.0 | 1.4 | $-396$
half-FM | 0.0 | 1.5 | $-381$
FiM | 2.1 | -1.4 | $-387$
AFM-V | 0.1 | $\pm 1.4$ | $-385$
Fe-cb |  $\pm$2.0 | 0.2 | $-219$
V-cb | 2.0 |  $\pm$1.2 | -237
V-cb | 0.1 |  $\pm$1.2 | -232
V-cb + Fe-stripes |  $\pm$2.2 | $\pm 1.2$ | $-409$

I thank W.E. Pickett and D. Parker for stimulating discussions related to this
work. I also acknowledge funding from the Office of Naval Research.

After this paper was accepted for publication, I became aware of another band
structure calculationChina . These author have considered the “Shein-
Ivanovskii” half-FM states and two antiferromagnetic states, with the
checkerboard (Neel) and stripe ordering in the Fe subslattice, and
unspecified, presumably ferromagnetic, ordering in the V subsystem. As clear
from the above, neither in this states represents an energy minimum even
within the corresponding symmetry group, therefore these authors arrived to an
incorrect conclusion that the lowest-energy magnetic state is characterized by
Neel order in the Fe subsystem.

## References

  * (1) X. Zhu, F. Han, G. Mu, P. Cheng, B. Shen, B. Zeng, and H.-H. Wen, Phys. Rev. B 79, 220512(R) (2009)
  * (2) I.I. Mazin and J. Schmalian, Physica C, 469, 614 (2009)
  * (3) V.I. Anisimov, E.Z. Kurmaev, A. Moewes, I.A. Izyumov, Physica C, 469, 442-447 (2009).
  * (4) W. L. Yang, A. P. Sorini, C-C. Chen, B. Moritz, W.-S. Lee, F. Vernay, P. Olalde-Velasco, J. D. Denlinger, B. Delley, J.-H. Chu, J. G. Analytis, I. R. Fisher, Z. A. Ren, J. Yang, W. Lu, Z. X. Zhao, J. van den Brink, Z. Hussain, Z.-X. Shen, and T. P. Devereaux. Phys. Rev. B80, 014508 (2009).
  * (5) K.-W. Lee and W. E. Pickett, arXiv:0908.2698 (unpublished)
  * (6) I used the standard LAPW code as implemented in the WIEN2k package with generalized gradient corrections. The orbitals weight used for the rest of the calculations are the LAPW orbital projection inside the corresponding muffin-tin spheres.
  * (7) M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. Parcollet, T. Miyake, A. Georges, and S. Biermann, Phys. Rev. B 80, 085101 (2009)
  * (8) Unfortunately, it is not clear at all whether the LDA+U treatment would be appropriate for V in this fluctuating system even if V itself is subject to strong Hubbard correlations.
  * (9) I. R. Shein and A. L. Ivanovskii, arXiv:0904.2671 (unpublished)
  * (10) K.-W. Lee and W. E. Pickett, private communication.
  * (11) G. Wang, M. Zhang, L. Zheng, and Z. Yang, Phys. Rev. B80, 184501 (2009).