# Probing many-body localization by excited-state VQE

Shuo Liu  Institute for Advanced Study, Tsinghua University, Beijing 100084,
China  Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong 518057, China
Shi-Xin Zhang<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua
University, Beijing 100084, China  Tencent Quantum Laboratory, Tencent,
Shenzhen, Guangdong 518057, China  Chang-Yu Hsieh<EMAIL_ADDRESS>Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong 518057, China
Shengyu Zhang  Tencent Quantum Laboratory, Tencent, Shenzhen, Guangdong
518057, China  Hong Yao<EMAIL_ADDRESS>Institute for Advanced Study,
Tsinghua University, Beijing 100084, China

###### Abstract

Non-equilibrium physics including many-body localization (MBL) has attracted
increasing attentions, but theoretical approaches of reliably studying non-
equilibrium properties remain quite limited. In this Letter, we propose a
systematic approach to probe MBL phases on a digital quantum computer via the
excited-state variational quantum eigensolver (VQE) and demonstrate convincing
results of MBL on a quantum hardware, which we believe paves a promising way
for future simulations of non-equilibrium systems beyond the reach of
classical computations in the noisy intermediate-scale quantum (NISQ) era.
Moreover, the MBL probing protocol based on excited-state VQE is NISQ-
friendly, as it can successfully differentiate the MBL phase from thermal
phases with relatively-shallow quantum circuits, and it is also robust against
the effect of quantum noises.

Introduction: Many-body localization (MBL) is a novel dynamical phenomenon
occuring in isolated many-body quantum systems. It has long been established
that the quantum systems may enter MBL phases in the presence of sufficiently
strong random disorder [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] or quasi-
periodic (QP) potential [13, 14, 15, 16, 17, 18] in one-dimensional (1D)
systems. In MBL phases, the system fails to thermally equilibrate and exhibits
exotic behaviors, such as “area law” entanglement for highly excited states
[9, 10, 11], logarithmic spread of entanglement [8], and emergent local
integrals of motion [11, 12]. Such exotic non-equilibrium phases are
qualitatively different from thermal phases that are often associated with the
eigenstate thermalization hypothesis (ETH) [19, 20, 21, 22, 23] and exhibit
“volume law” entanglement in highly excited states [24, 25]. However, various
aspects of MBL remain elusive so far; for instance, whether MBL phases can
survive in more than one dimensions is still under debate considering non-
perturbative avalanche mechanism [26, 27, 28].

To advance our understanding of non-equilibrium quantum phases, it is crucial
to unambiguously detect and characterize possible MBL phases for various
systems. Numerically, one can probe MBL by calculating the entanglement
entropy or level statistics for eigenstates of a Hamiltonian via exact
diagonalization, or directly simulating dynamical signatures such as charge
imbalance or logarithmic entanglement spreading, following a quantum quench
with time-evolution methods. But, existing numerical approaches are often
severely restricted due to the exponential scaling of Hilbert space as well as
the excessively long evolution time required to simulate steady-state
behaviors. Turning to experimental investigations, current hardware platforms
still face many challenges, such as short coherence time and limited
controllability for the Hamiltonian engineering. That whether MBL exists in
more general systems, for example, higher dimensional systems [29, 30, 31,
32], has not been rigorously established from current numerical simulations or
analog-simulation based experiments.

Quantum computer naturally helps as it can potentially simulate large quantum
systems, beyond the capability of classical computers. A common approach is to
directly simulate the time evolution of a quantum system using a quantum
circuit and determine whether the system will thermalize after a long time by
measuring dynamical observables such as charge imbalance. However, currently
such idea is not NISQ-friendly since the required simulation time to confirm
the existence of MBL can be much longer than the coherence time of the
currently available quantum hardware. Namely, we are faced with the same
challenge as other analog experiments without quantum error corrections.
Instead, here we propose a more NISQ-friendly approach to detect MBL, relying
less on quantum hardware resources and being more robust against noises in
quantum hardware.

In this Letter, we propose a general MBL-probing protocol basing on
variational quantum algorithms (VQAs). Recently, various quantum-classical
hybrid variational algorithms, such as variational quantum eigensolver [33]
and quantum approximate optimization algorithm (QAOA) [34, 35, 36, 37], which
are tailored for the NISQ hardware [38, 39, 40], have been proposed. VQE is
not only a representative VQA but also holds great potentials for near-term
applications. For different practical problems encoded in VQE, we can define
different objective functions (also known as cost functions). And the
solutions are expected to give the minimum objective function. For ground
state (lowest eigenstate) preparation for a Hamiltonian $H$, the objective
function is defined as the expectation value of $H$. Firstly, we will generate
a trial state $|\psi(\theta)\rangle=U(\theta)|\psi_{0}\rangle$ to approximate
the solution, where $|\psi_{0}\rangle$ is a given initial state and
$U(\theta)$ is an unitary matrix represented by the parameterized quantum
circuit (PQC) as shown in Fig. 1. When the PQC is deep enough, the ground
state can be exactly written as $|\psi(\theta)\rangle$ for suitable
parameters. Then, by updating the parameters in shallow PQC based on gradient
descent, we will find a final converged state $|\psi(\theta)\rangle$ whose
objective function can not decrease anymore and $|\psi(\theta)\rangle$ will be
the approximate ground state. The VQE has been exploited in a variety of
contexts from quantum chemistry [41, 42, 43, 44, 45, 46], many-body physics
[47, 48, 49, 35, 50, 51, 52, 53], to lattice gauge theories [54, 55].

VQE can also be customized to search for excited states. There are several
proposals utilizing VQE to discover low-lying excited states of quantum many-
body systems, such as orthogonality constrained VQE [56, 57, 58, 59], which
adds penalty projector terms to the Hamiltonian cost function that project out
lower energy states, subspace expansion method [60, 61], which prepares
variational states that span the low-energy manifold. Furthermore, VQE can be
adapted to search for highly excited states. The objective functions utilized
in this case are either the energy variance $C(\theta)=\langle
H^{2}\rangle-\langle H\rangle^{2}$ [62, 63, 64, 65, 66, 67, 68, 69, 70], which
only vanishes for eigenstates of $H$ and is greater than zero for the
superposition of different eigenstates of $H$ (see the Supplemental Material
for details), or $\langle(H-\lambda)^{2}\rangle$ [71, 72, 33, 73, 74], which
specifically targets an excited state closest to the energy $\lambda$. The
many-body spectrum can be reconstructed by scanning through a range of
$\lambda$ values within the energy width of $H$.

In this Letter, we propose to use the excited-state VQE with the energy-
variance objective function to detect MBL phases. If the circuit ansatz is
deep enough and with sufficient expressivity, then the converged state should
be an eigenstate of the system. One can then measure the quantum state to
infer purity or real-space inverse participation ratio (IPR) [75, 76] as the
indicators of MBL phases. This option is straightforward, but it requires a
(possibly exponentially) deep and coherent quantum circuit beyond the NISQ
regime. We propose an alternative approach, where only the excited-state VQE
with a NISQ compatible ansatz is required. In this case, the converged state
is not necessarily an eigenstate, especially on the thermal side, where the
eigenstate manifests “volume law” entanglement and may not be fully
represented with a shallow circuit. Therefore, the indicator to differentiate
between MBL and ergodic phases is replaced with the converged performance of
the excited-state VQE. On the MBL side, the final state is more quickly
converged in terms of the eigenspace distribution. Accordingly, we propose an
experimentally measurable quantity that could help assess a state’s
convergence in the eigenspace. Since our proposed method does not rely on any
properties of a particular model, it can be readily applied to a broad range
of systems.

Figure 1: The circuit structure for the excited-state VQE and eigenstate
witness measurement. In preparation circuit part, we first prepare an
antiferromagnetic (AF) state $|0101...\rangle$ from the initial state
$|0000...\rangle$ and then apply one layer of two-qubit entangling gates to
generate the input state $|\psi_{0}\rangle$ for the PQC. The layer of two-
qubit entangling gates is denoted as $U_{0}(\vec{\theta}^{0})$ and
$U_{0}(\vec{\theta}^{0})$ =
$\prod_{j=0}^{N-1}\mathrm{exp}(\frac{i\pi\theta^{0}_{j}}{4}(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}))$.
We assume two-qubit entangling gates on different qubits sharing the same
weight $\theta^{0}$ for simplicity. In the PQC part, each block of the ansatz
is defined as
$U(\vec{\theta}^{k})=(\prod_{j=0}^{N-1}\mathrm{Rz}(\theta^{k}_{2,j}))U_{0}(\vec{\theta}^{k}_{1})$,
which respects $U(1)$ symmetry. And two-qubit entangling gates on different
qubits have different weights. The depth of the PQC (number of blocks)
$\mathrm{D}_{\mathrm{VQE}}$ can be adjusted.

Model: We illustrate the proposed method by investigating the interacting
Aubry-André (AA) model [77, 78, 79, 14, 80, 16], a well-studied system hosting
the many-body localization transition. The Hamiltonian of the interacting AA
model reads:

$\displaystyle
H=\sum_{i}(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}+V_{0}\sigma_{i}^{z}\sigma_{i+1}^{z})$
$\displaystyle~{}~{}~{}~{}~{}~{}+W\sum_{i=1}^{N}\cos(2\pi\eta
i+\phi)\sigma_{i}^{z},$ (1)

where $W$ is the strength of quasi-periodic potential, $\sigma^{\alpha}$ are
Pauli matrices, $N$ is the size of the system, and $\phi$ is the phase of the
cosine potential. We set $V_{0}=0.5$ and $\eta=(\sqrt{5}-1)/2$ throughout the
work. We can numerically determine the MBL transition point via the level
spacing ratios for this model [3, 81] (see the Supplemental Material for
details).

Circuit ansatz: Our proposed circuit ansatz for the excited-state VQE consists
of two parts: the input-state preparation circuit and the parameterized
quantum circuit, acting as the variational ansatz to optimize the cost
function, as shown in Fig. 1. Since the AA model conserves the total spin
polarization along the $z$ direction, we focus on the total spin
$M_{z}=\sum_{i}\sigma_{i}^{z}=0$ sector. Namely, the quantum gates employed in
the circuit ansatz should respect this $U(1)$ symmetry. The cost function for
the excited-state VQE is the energy variance $\langle H^{2}\rangle-\langle
H\rangle^{2}$, as mentioned before.

Figure 2: EIPR of the input state $|\psi_{0}\rangle$ obtained from the
preparation circuit for $W=1.5$ and $8.0$, respectively. Here $N=12$.

Eigenspace inverse participation ratio: The final converged state
$|\psi\rangle$ obtained from the excited-state VQE is a superposition of
eigenstates, localized within an energy window, which may not be easily
resolved by the energy variance indicator because of the extremely small
energy separation between adjacent energy levels (see the Supplemental
Material for details). Theoretically, we may use the eigenspace inverse
participation ratio (EIPR) [74] to gauge the extent of convergence onto an
eigenstate under the excited-state VQE minimization:

$\displaystyle\mathrm{EIPR}(|\psi\rangle)=\sum_{n}|\psi_{n}|^{4},$ (2)

where $|\psi\rangle=\sum_{n}\psi_{n}|n\rangle$ and
$H|n\rangle=\lambda_{n}|n\rangle$ with $|n\rangle$ is the $n$-th eigenstate of
the Hamiltonian $H$. For each exact eigenstate of $H$, its EIPR is one. For a
state which is a linear superposition of different eigenstates, its EIPR is
less than one. For the maximally randomized state in eigenspace
$|\psi\rangle=\sum_{n}\frac{1}{\sqrt{2^{N}}}|n\rangle$, where $2^{N}$ is the
dimension of the Hilbert space of the system;
$\mathrm{EIPR}=1/2^{N}\rightarrow 0$ as $N\rightarrow\infty$ in the
thermodynamic limit. Consequently, EIPR is capable in efficiently determining
whether the final converged state is exactly an eigenstate of the Hamiltonian
$H$ or a superposition of multiple eigenstates.

We firstly investigate the EIPR for the input state $|\psi_{0}\rangle$,
created with the preparation circuit of Fig. 1. The results are shown in Fig.
2. When the system is in the MBL phase ($W=8.0$), some of the shallow-circuit
input states (adjusted by the parameter $\theta^{0}$) may achieve high EIPR,
implying the “area law” entanglement [8, 9, 10, 11] for highly-excited
eigenstates. On the other hand, when the system is in the thermal phase
($W=1.5$), all input states consistently score low EIPR, indicating that a
shallow circuit has difficulty to approximate the system’s eigenstates.
Overall, the results shown in Fig. 2 are consistent with the understanding of
the AA model. Furthermore, we find a positive correlation between EIPR of
input states and EIPR of output states, as well as a positive correlation
between energy of input states and energy of output states (see the
Supplemental Material for details). For the following investigations, we then
fix $\theta^{0}=0.4$ for input-state preparation. As shown in Fig. 2, the MBL
phase has a smaller initial EIPR and the input state’s energy is in the middle
of energy spectrum in both phases. By this choice, under the excited-state
VQE, the output states will likely to converge to some highly excited states,
and a high EIPR for output state in MBL phase can be reliably attributed to
the nature of the phase instead of a purely better optimization start point.

With the noiseless quantum circuit simulation [82], we calculate EIPR for
converged states with varying PQC depth. Theoretically, the expressivity of
the variational ansatz increases with the circuit depth, and the achieved EIPR
should also improve correspondingly. Indeed, the results for converged states’
EIPR, as shown in Fig. 3, support this intuition. Although the initial EIPR of
the MBL phase is smaller than that of the thermal phase and the optimization
starts from the same input state, the final EIPR of the MBL phase is much
larger and extremely close to $1$ even with a relatively shallow quantum
circuit. There is a clear gap of EIPR between thermal and MBL phases, and the
gap grows with system size (see the Supplemental Material for details). These
results confirm that excited-state VQE with shallow circuits indeed performs
qualitatively better for MBL systems.

Figure 3: EIPR of final converged states from excited-state VQE on 12-qubit
system with different potential strength $W$. Figure 4: (a),
$\ln(1-\mathrm{EIPR})$ of final converged states. (b), $\ln(1-r)$ of final
converged states, where $r$ is the eigenstate witness and experimentally
accessible. Here $N=12$.

Experimental relevance: Though EIPR of the output states significantly differs
between MBL and thermal phases in numerical simulations, it can not be
directly observed in experiments. To determine whether a state is sufficiently
enough to an eigenstate, an experimentally accessible proxy is desired. To
this end, the eigenstate witness [74] is such a quantity which can be
experimentally determined via purity measurements. The circuit for eigenstate
witness measurement is shown in Fig. 1.

We illustrate the notation of the eigenstate witness as follows. The output
state obtained from the excited-state VQE can be written as
$|\psi\rangle=\sum_{n}\alpha_{n}|n\rangle$, where $|n\rangle$ represents
eigenstate of $H$ with eigenenergy $\lambda_{n}$. With an extra ancilla qubit
initialized in a superposition state
$|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$, after a controlled time evolution
under $H$, the reduced density matrix for the ancilla qubit reads:

$\displaystyle\rho_{\mathrm{reduced}}=\left(\begin{matrix}\frac{1}{2}&\frac{1}{2}\sum_{n}|\alpha_{n}|^{2}\mathrm{e}^{i\lambda_{n}t}\\\
\frac{1}{2}\sum_{n}|\alpha_{n}|^{2}\mathrm{e}^{-i\lambda_{n}t}&\frac{1}{2}\end{matrix}\right),$
(3)

and the eigenstate witness is defined as the purity of the ancilla qubit:
$r=\mathrm{Tr}(\rho_{\mathrm{reduced}}^{2}),$ which can be estimated using
randomized measurements in experiments [83, 84, 85] or a simple state
tomography. Such witness is lower bounded by the corresponding EIPR of the
output states and can reflect the convergence performance of the excited-state
VQE (see the Supplemental Material for details).

We choose the evolution time $t=1.0/W$, as the bandwidth is roughly
proportional to $W$. The results of $\mathrm{ln}(1-r)$ are shown in Fig. 4.
For each data point in Fig. 4, the number is averaged over 10 best out of 100
independent optimization results for the excited-state VQE. Both EIPR and $r$
of the MBL phase are larger than those of the thermal phase. These differences
will be more prominent as the system size increases.

Theoretically, $\ln(1-r)$ should approach the minus infinity if the output
state is exactly an eigenstate of the system. However, due to the
effectiveness of an optimization routine and numerical accuracy of the
simulation, the $(1-r)$ value may only reach an “effective zero” from above.
In Fig. 4, $(1-\mathrm{EIPR})$ and $(1-r)$ have this “effective zero” which is
roughly determined by the data from $W=100.0$, supposedly in a deep region of
the MBL phase. The effective zero of $(1-r)$ is around $-8.9$ in the log
scale. And the circuit depth required to reach this effective zero depends on
Hamiltonian parameter $W$. We estimate the required depth in Table. 1 (the
required depth is fitted by assuming linear relation between $\ln(1-r)$ and
the circuit depth). On the thermal side, the depth required for perfect
excited-state convergence is close to $2^{N}$; while on the MBL side, the
depth required is in the order of $O(N)$. The depths needed to reach effective
zero can serve as an another indicator differentiating between the MBL and
thermal phases in systems whose phases are unknown. In particular, it is
easier to distinguish the two phases when the system size increases.

W |  1.5 |  2.0 |  2.5 |  3.5 |  4.0 |  8.0
---|---|---|---|---|---|---
Depths | 245 | 197 | 256 | 52 | 38 | 14
Table 1: Fitted circuit depths required to perfectly converge to the excited
state for different $W$. Figure 5: (a), $r$ of final converged states from a
noisy simulation without Trotter decomposition. (b), $r$ of final converged
states from a noisy simulation with one time slice Trotter decomposition. Here
$N=8$.

Effect of quantum noise: Our approach of probing MBL phases has potential
quantum advantage as it is scalable to larger systems and higher dimensions.
However, for current quantum computing hardware, we must consider the effects
of noise which may compromise the performance of excited-state VQE. One
important question is whether we can apply our method to a real device with
noises. To answer this question, we carry out similar computation in the
presence of a quantum depolarizing channel with noise strength $p=10^{-3}$
after each two-qubit gate. The results are averaged over $5$ best of $20$
independent trials on an 8-qubit system (the results for larger size systems
with the presence of quantum noise obtained by the Monte Carlo trajectory
method can be found in the Supplemental Material). As indicated in Fig. 5(a),
$r$ does not keep increasing with PQC depths any more; instead, it could
decrease due to the accumulated noisy effects for deeper PQC. Nonetheless, the
eigenstate witness $r$ behaves sufficiently distinct in the two phases. The
magnitude of $r$ from the MBL side is still significantly larger than that
from the thermal side. Besides, $r$ shows better noise resilience when the
system enters the MBL side. In the MBL phase, local perturbations spread only
logarithmically in time [86], as opposed to the algebraic spreading in
thermalizing dynamics. Therefore, $r$ of the thermal phase gets more severely
affected due to the faster spreading of quantum noise.

Figure 6: Eigenstate witness measured on the noisy simulator and real hardware
(IBM_Santiago) (in the inset ∗ means the results are processed with readout
error mitigation). The estimation error bars due to finite shots of
measurement are also included.

When implementing the controlled time evolution in a real device, one also
needs to consider the error brought by the Trotter decomposition. Here we
design a type of decomposition that reduces the number of two-qubit gates, and
we propagate the system with only one time slice to further reduce the total
circuit depth for the digital quantum simulation (see the Supplemental
Material for details). The eigenstate witness $r$, after incorporating the
effects of both noise and Trotter decomposition, is shown in Fig. 5(b), which
still behaves distinctly in the two phases. The numerical results confirm that
our proposal is NISQ-friendly and ready to be validated in a quantum device.

Real Hardware Experiments: We apply our method to a four-site model in the
available quantum hardware. When the quantum computer is ready for simulating
large enough systems, the benefits of the current proposal will become more
evident. Since the controlled time evolution module requires a large number of
two-qubit gates, we further reduce the depth by utilizing a variational
quantum circuit to approximate the controlled time evolution module. Such a
variational circuit is based on the hardware-efficient ansatz and has less
number of two-qubit gates in total while maintaining a high fidelity for the
exact controlled time evolution (see the Supplemental Material for more detail
of real hardware experiments).

We demonstrate our method in a noisy hardware simulator based on IBMQ hardware
device (specifically we use the IBM_Santiago instance). The results are shown
in Fig. 6. A qualitative difference is clearly seen which is sufficient to
distinguish the two phases. In fact, for this case, the difference obtained in
quantum hardware experiments is even more evident.

Discussions and concluding remarks: We have demonstrated that the excited-
state VQE with EIPR constitutes a reliable method to probe MBL phase. Since
this proposed method requires only a shallow circuit, it can be readily
executed on a NISQ hardware when we substitute the EIPR with the eigenstate
witness protocol. Compared to the conventional approaches that characterize
the MBL phase with late-time dynamics, the evolution time required in our
method is much shorter and compatible with the current-generation quantum
computers.

There are many promising directions for further studies. In the present work
we have considered only one circuit ansatz. Other different PQC ansatz for
excited-state VQE can also be systematically examined for potential benefits.
The circuit ansatz with better expressiveness may characterize the highly
excited states better, and the optimal ansatz can be automatically designed
via techniques from quantum architecture search [87, 88]. We can also utilize
post-processing enhancements [89] which hopefully improve the performance of
excited-state VQE. As our method is independent of the microscopic details of
the model considered in this work, it can be straightforwardly extended to
investigate other models which might host many-body localization transitions.
Moreover, our method can hopefully investigate whether MBL phases can survive
in higher dimensions in the future, once the quantum hardware can offer longer
coherence time.

Acknowledgements: We thank Zhou-Quan Wan for helpful discussions. This work is
supported in part by the NSFC under Grant No. 11825404 (SXZ, SL, and HY), the
CAS Strategic Priority Research Program under Grant No. XDB28000000 (HY), and
Beijing Municipal Science and Technology Commission under Grant No.
Z181100004218001 (HY).

## References

  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006).
  * Oganesyan and Huse [2007a] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Physical Review B 75, 155111 (2007a).
  * Oganesyan and Huse [2007b] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007b).
  * Znidaric  _et al._ [2008] M. Znidaric, T. c. v. Prosen, and P. Prelovsek, Many-body localization in the heisenberg $xxz$ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Monthus and Garel [2010] C. Monthus and T. Garel, Many-body localization transition in a lattice model of interacting fermions: Statistics of renormalized hoppings in configuration space, Phys. Rev. B 81, 134202 (2010).
  * Cuevas  _et al._ [2012] E. Cuevas, M. Feigel’man, L. Ioffe, and M. Mezard, Level statistics of disordered spin-1/2 systems and its implications for materials with localized Cooper pairs, Nature Communications 3, 1128 (2012).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Vosk and Altman [2013] R. Vosk and E. Altman, Many-body localization in one dimension as a dynamical renormalization group fixed point, Phys. Rev. Lett. 110, 067204 (2013).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013a).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013b).
  * Huse and Oganesyan [2014] D. A. Huse and V. Oganesyan, A phenomenology of certain many-body-localized systems, Physical Review B 90, 174202 (2014).
  * Iyer  _et al._ [2013] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Many-body localization in a quasiperiodic system, Phys. Rev. B 87, 134202 (2013).
  * Modak and Mukerjee [2015] R. Modak and S. Mukerjee, Many body localization in the presence of a single particle mobility edge, Physical Review Letters 115, 230401 (2015).
  * Zhang and Yao [2018] S.-X. Zhang and H. Yao, Universal properties of many-body localization transitions in quasiperiodic systems, Physical Review Letters 121, 206601 (2018).
  * Kohlert  _et al._ [2019] T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. D. Sarma, I. Bloch, and M. Aidelsburger, Observation of many-body localization in a one-dimensional system with single-particle mobility edge, arXiv:1809.04055 (2019).
  * Zhang and Yao [2019] S.-X. Zhang and H. Yao, Strong and Weak Many-Body Localizations, arXiv:1906.00971 (2019).
  * Ghosh  _et al._ [2020] S. Ghosh, J. Gidugu, and S. Mukerjee, Transport in the non-ergodic extended phase of interacting quasiperiodic systems, Physical Review B 102, 224203 (2020).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From Quantum Chaos and Eigenstate Thermalization to Statistical Mechanics and Thermodynamics, Advances in Physics 65, 239 (2016).
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Reports on Progress in Physics 81, 082001 (2018).
  * Deutsch [2010] J. M. Deutsch, Thermodynamic entropy of a many-body energy eigenstate, New Journal of Physics 12, 075021 (2010).
  * Garrison and Grover [2018] J. R. Garrison and T. Grover, Does a single eigenstate encode the full hamiltonian?, Phys. Rev. X 8, 021026 (2018).
  * De Roeck and Huveneers [2017] W. De Roeck and F. m. c. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Thiery  _et al._ [2018] T. Thiery, F. m. c. Huveneers, M. Müller, and W. De Roeck, Many-body delocalization as a quantum avalanche, Phys. Rev. Lett. 121, 140601 (2018).
  * Potirniche  _et al._ [2019] I.-D. Potirniche, S. Banerjee, and E. Altman, Exploration of the stability of many-body localization in $d>1$, Phys. Rev. B 99, 205149 (2019).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Bordia  _et al._ [2016] P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider, Coupling Identical one-dimensional Many-Body Localized Systems, Physical Review Letters 116, 140401 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Probing Slow Relaxation and Many-Body Localization in Two-Dimensional Quasiperiodic Systems, Physical Review X 7, 041047 (2017).
  * Wahl  _et al._ [2019] T. B. Wahl, A. Pal, and S. H. Simon, Signatures of the Many-body Localized Regime in Two Dimensions, Nature Physics 15, 164 (2019).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a quantum processor, Nature Communications 5, 4213 (2014).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv:1411.4028 (2014).
  * McClean  _et al._ [2016a] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016a).
  * Wang  _et al._ [2018] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, Quantum Approximate Optimization Algorithm for MaxCut: A Fermionic View, Physical Review A 97, 022304 (2018).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Physical Review X 10, 021067 (2020).
  * Cerezo  _et al._ [2020] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational Quantum Algorithms, arXiv:2012.09265 (2020).
  * Bharti  _et al._ [2021] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum (NISQ) algorithms, arXiv:2101.08448 (2021).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, Journal of the Physical Society of Japan 90, 032001 (2021).
  * 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, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, A. Megrant, J. Y. Mutus, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Scalable Quantum Simulation of Molecular Energies, Physical Review X 6, 031007 (2016).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient Variational Quantum Eigensolver for Small Molecules and Quantum Magnets, Nature 549, 242 (2017).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Computation of molecular spectra on a quantum processor with an error-resilient algorithm, Phys. Rev. X 8, 011021 (2018).
  * Armaos  _et al._ [2019] V. Armaos, D. A. Badounas, and P. Deligiannis, Computational Chemistry on Quantum Computers: Ground state estimation, arXiv:1907.00362 (2019).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019).
  * Tilly  _et al._ [2020] J. Tilly, G. Jones, H. Chen, L. Wossnig, and E. Grant, Computation of molecular excited states on IBM quantum computers using a discriminative variational quantum eigensolver, Physical Review A 102, 062425 (2020).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Phys. Rev. A 92, 042303 (2015).
  * Dallaire-Demers  _et al._ [2018] P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, Low-depth circuit ansatz for preparing correlated fermionic states on a quantum computer, arXiv:1801.01053 (2018).
  * Ho and Hsieh [2019] W. W. Ho and T. H. Hsieh, Efficient variational simulation of non-trivial quantum states, SciPost Phys. 6, 29 (2019).
  * Liu  _et al._ [2019] J.-G. Liu, Y.-H. Zhang, Y. Wan, and L. Wang, Variational Quantum Eigensolver with Fewer Qubits, Physical Review Research 1, 023025 (2019).
  * Cai [2020] Z. Cai, Resource Estimation for Quantum Variational Simulations of the Hubbard Model, Physical Review Applied 14, 014059 (2020).
  * Uvarov  _et al._ [2020] A. Uvarov, J. Biamonte, and D. Yudin, Variational Quantum Eigensolver for Frustrated Quantum Systems, Physical Review B 102, 075104 (2020).
  * Tamiya  _et al._ [2021] S. Tamiya, S. Koh, and Y. O. Nakagawa, Calculating nonadiabatic couplings and Berry’s phase by variational quantum eigensolvers, Physical Review Research 3, 023244 (2021).
  * Klco  _et al._ [2018] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, and M. J. Savage, Quantum-classical computation of schwinger model dynamics using quantum computers, Phys. Rev. A 98, 032331 (2018).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller, Self-Verifying Variational Quantum Simulation of the Lattice Schwinger Model, Nature 569, 355 (2019).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Quantum fingerprinting, Phys. Rev. Lett. 87, 167902 (2001).
  * Cincio  _et al._ [2018] L. Cincio, Y. Subası, A. T. Sornborger, and P. J. Coles, Learning the quantum algorithm for state overlap, New Journal of Physics 20, 113022 (2018).
  * Higgott  _et al._ [2019] O. Higgott, D. Wang, and S. Brierley, Variational Quantum Computation of Excited States, Quantum 3, 156 (2019).
  * Jones  _et al._ [2019] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, Variational quantum algorithms for discovering hamiltonian spectra, Phys. Rev. A 99, 062304 (2019).
  * McClean  _et al._ [2017] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states, Phys. Rev. A 95, 042308 (2017).
  * Nakanishi  _et al._ [2019] K. M. Nakanishi, K. Mitarai, and K. Fujii, Subspace-search variational quantum eigensolver for excited states, Physical Review Research 1, 033062 (2019).
  * Bartlett  _et al._ [1935] J. H. Bartlett, J. J. Gibbons, and C. G. Dunn, The Normal Helium Atom, Physical Review 47, 679 (1935).
  * Umrigar  _et al._ [1988] C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Optimized trial wave functions for quantum Monte Carlo calculations, Physical Review Letters 60, 1719 (1988).
  * Marotta and Siringo [2005] L. Marotta and F. Siringo, A variational method from the variance of energy, The European Physical Journal C 44, 293 (2005).
  * Umrigar and Filippi [2005] C. J. Umrigar and C. Filippi, Energy and Variance Optimization of Many-Body Wave Functions, Physical Review Letters 94, 150201 (2005).
  * Khemani  _et al._ [2016] V. Khemani, F. Pollmann, and S. Sondhi, Obtaining Highly Excited Eigenstates of Many-Body Localized Hamiltonians by the Density Matrix Renormalization Group Approach, Physical Review Letters 116, 247204 (2016).
  * Pollmann  _et al._ [2016] F. Pollmann, V. Khemani, J. I. Cirac, and S. L. Sondhi, Efficient variational diagonalization of fully many-body localized Hamiltonians, Physical Review B 94, 041116 (2016).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational Neural-Network Ansatz for Steady States in Open Quantum Systems, Physical Review Letters 122, 250503 (2019).
  * Zhang  _et al._ [2020a] D.-B. Zhang, Z.-H. Yuan, and T. Yin, Variational quantum eigensolvers by variance minimization, arXiv:2006.15781 (2020a).
  * Zhang  _et al._ [2021a] F. Zhang, N. Gomes, Y. Yao, P. P. Orth, and T. Iadecola, Adaptive variational quantum eigensolvers for highly excited states, arXiv:2104.12636 (2021a).
  * MacDonald [1934] J. K. L. MacDonald, On the Modified Ritz Variation Method, Physical Review 46, 828 (1934).
  * Wang and Zunger [1994] L. Wang and A. Zunger, Solving Schrödinger’s equation around a desired energy: Application to silicon quantum dots, The Journal of Chemical Physics 100, 2394 (1994).
  * McClean  _et al._ [2016b] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016b).
  * Santagati  _et al._ [2018] R. Santagati, J. Wang, A. A. Gentile, S. Paesani, N. Wiebe, J. R. McClean, S. Morley-Short, P. J. Shadbolt, D. Bonneau, J. W. Silverstone, D. P. Tew, X. Zhou, J. L. O’Brien, and M. G. Thompson, Witnessing eigenstates for quantum simulation of Hamiltonian spectra, Science Advances 4 (2018).
  * Plerou  _et al._ [1999] V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Amaral, and H. E. Stanley, Universal and non-universal properties of cross-correlations in financial time series, Physical Review Letters 83, 1471 (1999).
  * Pradhan and Sridhar [2000] P. Pradhan and S. Sridhar, Correlations due to localization in quantum eigenfunctions of disordered microwave cavities, Physical Review Letters 85, 2360 (2000).
  * [77] S. Aubry and G. André, Analyticity breaking and anderson localization in incommensurate lattices, Ann. Israel Phys. Soc 3, 18.
  * Lahini  _et al._ [2009] Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, N. Davidson, and Y. Silberberg, Direct observation of a localization transition in quasi-periodic photonic lattices, Physical Review Letters 103, 013901 (2009).
  * Biddle  _et al._ [2009] J. Biddle, B. Wang, D. J. Priour Jr., and S. D. Sarma, Localization in one-dimensional incommensurate lattices beyond the Aubry-Andr\’e model, Physical Review A 80, 021603 (2009).
  * Dutta  _et al._ [2018] A. Dutta, S. Mukerjee, and K. Sengupta, Many-body localized phase of bosonic dipoles in a tilted optical lattice, Physical Review B 98, 144205 (2018).
  * Throckmorton and Sarma [2021] R. E. Throckmorton and S. D. Sarma, Studying many-body localization in exchange-coupled electron spin qubits using spin-spin correlations, Physical Review B 103, 165431 (2021).
  * Zhang  _et al._ [2022] S.-X. Zhang, J. Allcock, Z.-Q. Wan, S. Liu, J. Sun, H. Yu, X.-H. Yang, J. Qiu, Z. Ye, Y.-Q. Chen, C.-K. Lee, Y.-C. Zheng, S.-K. Jian, H. Yao, C.-Y. Hsieh, and S. Zhang, TensorCircuit: a Quantum Software Framework for the NISQ Era, arXiv:2205.10091 (2022).
  * Elben  _et al._ [2019] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A 99, 052323 (2019).
  * Elben  _et al._ [2020] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Cross-platform verification of intermediate scale quantum devices, Phys. Rev. Lett. 124, 010504 (2020).
  * Haug  _et al._ [2021] T. Haug, C. N. Self, and M. S. Kim, Large-scale quantum machine learning, arXiv:2108.01039 (2021).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Zhang  _et al._ [2020b] S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Differentiable Quantum Architecture Search, arXiv:2010.08561 (2020b).
  * Zhang  _et al._ [2021b] S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Neural Predictor based Quantum Architecture Search, arXiv:2103.06524 (2021b).
  * Zhang  _et al._ [2021c] S.-X. Zhang, Z.-Q. Wan, C.-K. Lee, C.-Y. Hsieh, S. Zhang, and H. Yao, Variational Quantum-Neural Hybrid Eigensolver, arXiv:2106.05105 (2021c).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Vatan and Williams [2004] F. Vatan and C. Williams, Optimal Quantum Circuits for General Two-Qubit Gates, Physical Review A 69, 032315 (2004).
  * Mandviwalla  _et al._ [2018] A. Mandviwalla, K. Ohshiro, and B. Ji, in 2018 IEEE International Conference on Big Data (Big Data) (IEEE, Seattle, WA, USA, 2018) pp. 2531–2537.

## Supplemental Materials

###  .1 A. Level spacing ratio for many-body localization transition

The level spacing ratio is defined as:

$\displaystyle
r_{\phi}=\frac{\mathrm{min}\\{\Delta^{n}_{\phi},\Delta^{n+1}_{\phi}\\}}{\mathrm{max}\\{\Delta^{n}_{\phi},\Delta^{n+1}_{\phi}\\}},$
(S1)

where $\Delta^{n}_{\phi}=E^{n}_{\phi}-E^{n+1}_{\phi}$ is the gap between
eigenenergy levels $n$ and $n+1$, and $\phi$ labels a given quasi-periodic
potential, $[r^{(n)}_{\phi}]$ is the level average of this ratio. In the
delocalized phase (thermal phase), the level spacings follow a Gaussian
orthogonal ensemble (GOE) distribution and the level-averaged level spacing
ratio is $[r^{(n)}_{\phi}]\approx 0.53$ [90]. In contrast, in the localized
phase, the level spacing follows a Poisson distribution, which gives
$[r^{(n)}_{\phi}]=2\ln 2-1\approx 0.39$ [90]. The results of level spacing
ratios for the interacting AA model are given in Fig. S1.

Figure S1: Exact diagonalization calculation: level spacing ratio with varying
quasi-periodic potential strength $W$. Each point is an average of 1000
different random phases $\phi$.

###  .2 B. The failure of energy variance to determine the VQE convergence

Consider a state $|\psi\rangle$:

$\displaystyle|\psi\rangle=\sum_{n}\psi_{n}|n\rangle,$ (S2)

where $|n\rangle$ is the $n$-th eigenstate of the Hamiltonian $H$ and
$H|n\rangle=\lambda_{n}|n\rangle$. The energy variance is

$\displaystyle\langle\psi|H^{2}|\psi\rangle-\langle\psi|H|\psi\rangle^{2}$
$\displaystyle=$ $\displaystyle\sum_{m,n}\psi^{*}_{m}\psi_{n}\langle
m|H^{2}|n\rangle-(\sum_{m,n}\psi^{*}_{m}\psi_{n}\langle m|H|n\rangle)^{2}$
$\displaystyle=$
$\displaystyle\sum_{n}|\psi_{n}|^{2}\lambda_{n}^{2}-(\sum_{n}|\psi_{n}|^{2}\lambda_{n})^{2}$
$\displaystyle=$
$\displaystyle(\sum_{m}|\psi_{m}|^{2})(\sum_{n}|\psi_{n}|^{2}\lambda_{n}^{2})-(\sum_{n}|\psi_{n}|^{2}\lambda_{n})^{2}$
$\displaystyle\geq$
$\displaystyle(\sum_{n}|\psi_{n}|^{2}\lambda_{n})^{2}-(\sum_{n}|\psi_{n}|^{2}\lambda_{n})^{2}$
$\displaystyle=$ $\displaystyle 0,$

where we have utilized the Cauchy–Schwarz inequality. When $|\psi\rangle$ is
an eigenstate of $H$, i.e., $|\psi\rangle=|n\rangle$, the energy variance is
zero. Otherwise, the energy variance is positive. Therefore, the energy
variance can be used as the cost function for the excited-state VQE.

The final converged state $|\psi\rangle$ obtained from the excited-state VQE
is a superposition of eigenstates, mainly localized within an energy window,
which may not be easily resolved by the energy variance indicator because of
the extremely small energy separation between adjacent energy levels. Assume
the final converged state is simply a superposition of two eigenstates:
$|\psi\rangle=a|E\rangle+b|E+\Delta\rangle$, where $a^{2}+b^{2}=1$.
$|E\rangle$ and $|E+\Delta\rangle$ are eigenstates with energy $E$ and
$E+\Delta$, respectively, and $\Delta$ is the energy gap which is extremely
close to zero. Then

$\displaystyle\langle H\rangle=\langle\psi|H|\psi\rangle$ $\displaystyle=$
$\displaystyle a^{2}E+b^{2}(E+\Delta)$ $\displaystyle=$ $\displaystyle
E+(1-a^{2})\Delta,$

and

$\displaystyle\langle H^{2}\rangle=\langle\psi|H^{2}|\psi\rangle$
$\displaystyle=$ $\displaystyle a^{2}E^{2}+b^{2}(E+\Delta)^{2}$
$\displaystyle=$ $\displaystyle E^{2}+(1-a^{2})(2E\Delta+\Delta^{2}).$

The energy variance is:

$\displaystyle\langle H^{2}\rangle-\langle
H\rangle^{2}=(a^{2}-a^{4})\Delta^{2}.$ (S6)

It will approach zero even if the final state is a superposition of different
eigenstates, as long as $\Delta$ is close to zero. In other words, the energy
variance is no longer a good indicator to verify whether the excited-state VQE
converged to one eigenstate in the limit of small $\Delta$.

###  .3 C. The EIPR and energy correlation between input states and output
states of VQE

The initial input state $|\psi_{0}\rangle$ is a superposition of eigenstates.
If the initial EIPR is large, which means one eigenstate has already been
dominant since the beginning, then the final converged state for a VQE
calculation shall flow to that eigenstate. As shown in Fig. S2, EIPR is
maximum when $\theta^{0}=1.0$ and $W=8.0$. We show the spectral decomposition
of this initial input state in the eigenspace in Fig. S3. The EIPR of the
final output state is almost $1$ and successfully converges to that dominant
eigenstate as shown in Fig. S3. Even if we randomize the initial parameters of
PQC in a small range, then all final converged states will still flow to that
eigenstate as shown in Fig. S4. When the initial EIPR is small, the input
state is more uniformly distributed in the eigenspace as shown in Fig. S5. In
such cases, the performance of the excited-state VQE deteriorates, the final
converged state is stuck as a superposition of different eigenstates.
Furthermore, the dominant eigenstate of the output state will change in
different independent trials as shown in Fig. S6.

Figure S2: EIPR of input state $|\psi_{0}\rangle$ determined by the
preparation circuit with $W=1.5,8.0$; $N=12$. Figure S3: Overlap with
different eigenstates: input state vs. final converged state: $W=8.0$; $N=10$;
depth=3. Figure S4: Final converged energy for $W=8.0$ and input state of good
$\theta^{0}$; $N=10$; depth=3. Figure S5: Overlap with different eigenstates:
input state vs. final converged state: $W=8.0$; $N=10$; depth=3. Figure S6:
Final converged energy for $W=8.0$ and input state of bad $\theta^{0}$;
$N=10$; depth=3.

Another interesting observation is that EIPR for the input states
$|\psi_{0}\rangle$, created with the preparation circuit, will decrease as the
system size increases as shown in Fig. S7. For a larger system, EIPR of
thermal phases and MBL phases are both close to $0$. Since our setup in this
paper starts from input state with small EIPR on purpose, the approach is
still suitable for larger size systems.

A straightforward question is whether the final converged states of the
thermal phase are low-lying excited states that also exhibit “area law”
entanglement, which may compromise the characteristic gap of EIPR between MBL
and thermal phases as shown in Fig. 3. As mentioned in the adaptive VQE-X
algorithm [70], we also find a weakly linear relationship between energy of
input states and energy of output states. With $\theta^{0}=0.4$, most of the
output states are highly excited states that exhibit “volume law” entanglement
(see Fig. S8). So the entanglement difference induces the gap of EIPR between
MBL and thermal phases.

Figure S7: EIPR of input state $|\psi_{0}\rangle$ determined by the
preparation circuit with $W=8.0$; $N=10,12,14$. Figure S8: Final output
states’ energy with $W=2.5$; $N=12$; depth=12. Most of the output states are
highly excited states.

###  .4 D. Scaling analysis of EIPR and eigenstate witness

The results of EIPR with $N=12$ in Fig. 3 show an identifiable gap between MBL
and thermal phases. To extend our method to larger size systems, we should
study the finite size scaling of this gap. To this end, we define:
$\Delta\mathrm{EIPR}=\mathrm{EIPR_{MBL}}-\mathrm{EIPR_{thermal}}$, where
$W=4.5$ for the MBL phase and $W=2.5$ for the thermal phase. As the system
size increases, we set PQC depth=$N$ that also scales with the system size. As
shown in Fig. S9, the gap of EIPR between MBL and thermal phases will be more
prominent for larger systems. As shown in the main text, on the MBL side, the
optimization depth required for the perfect excited-state convergence is in
the order $O(N)$. On the thermal side, the depth required is in the order of
$O(2^{N})$. The ideal gap is close to 1 when the optimization for thermal
phases completely fails with only depth=$N$ PQC for larger system.

For current quantum computing hardware, we must consider the effects of the
noise. In particular, the noise accumulates with the system size/circuit
depth. One important question is whether we can apply our method to a real
device with moderate noise. To answer this question, we have also studied the
finite size scaling of $\Delta r$: the eigenstate witness difference between
the MBL phase ($W=4.5$) and the thermal phase ($W=2.5$). To extend our method
to larger size systems, We carry out Monte Carlo trajectory simulation [82] in
the presence of a quantum depolarizing channel with noise strength $p=10^{-3}$
after each two-qubit gate. Quantum channels are added to the circuit by
providing a list of their associated Kraus operators, and the mixed state is
effectively simulated by an ensemble of pure states. Via this numerical
technique, larger size systems are accessible compared to density matrix
simulator. As the system size increases, we set PQC depth=$N/2$ that also
scales with the system size. As shown in Fig. S10, the eigenstate witness $r$
decreases with increasing system size due to stronger effective noise.
However, the eigenstate witness difference $\Delta r$ becomes larger when the
system size increases, remaining qualitatively the same as the noiseless case
shown in Fig. S9.

Based on these convincing numerical results, it is reliable to apply our
method to larger size systems on real noisy devices.

Figure S9: $\Delta\mathrm{EIPR}$ with $N=8,10,12,14$. Figure S10: (a), $r$ of
final converged states without Trotter decomposition. (b), $r$ of final
converged states with one time slice Trotter decomposition. The insets are the
eigenstate witness difference $\Delta r$ between the MBL phase and the thermal
phase. $N=6,8,10,12$.

###  .5 E. Eigenstate witness

The state after the controlled time evolution can be written as:

$\displaystyle\frac{|0\rangle|\psi\rangle+|1\rangle\mathrm{e}^{-iHt}|\psi\rangle}{\sqrt{2}}.$
(S7)

After a partial trace over the system qubit register, the reduced density
matrix for the ancilla qubit reads:

$\displaystyle\rho_{\mathrm{reduced}}=\left(\begin{matrix}\frac{1}{2}&\frac{1}{2}\sum_{n}|\alpha_{n}|^{2}\mathrm{e}^{i\lambda_{n}t}\\\
\frac{1}{2}\sum_{n}|\alpha_{n}|^{2}\mathrm{e}^{-i\lambda_{n}t}&\frac{1}{2}\end{matrix}\right),$
(S8)

where $|\psi\rangle=\sum_{n}\alpha_{n}|n\rangle$ such that
$H|n\rangle=\lambda_{n}|n\rangle$. The eigenstate witness is defined as:

$\displaystyle r=\mathrm{Tr}(\rho_{\mathrm{reduced}}^{2})$
$\displaystyle=\frac{1}{2}(\sum_{n}|\alpha_{n}|^{4})$ (S9)
$\displaystyle+\frac{1}{2}(1+\sum_{n\neq
m}|\alpha_{n}|^{2}|\alpha_{m}|^{2}\cos((\lambda_{n}-\lambda_{m})t))$
$\displaystyle\geq\sum_{n}|\alpha_{n}|^{4}$ $\displaystyle=\mathrm{EIPR}.$

We can see the eigenstate witness $r$ is lower bounded by EIPR. For almost all
$t$, obtaining an eigenstate witness of $1$ implies that $|\psi\rangle$ is an
eigenstate of the Hamiltonian $H$. When the evolution time is assumed to be
long, $r=\frac{1}{2}(1+\text{EIPR})$; when the evolution time is assumed to be
short, i.e. $\text{max}_{m,n}|\lambda_{m}-\lambda_{n}|t\leq\frac{\pi}{2}$,
then

$\displaystyle
r\leq\frac{1+\cos(\sum_{m,n}|\alpha_{m}|^{2}|\alpha_{n}|^{2}(\lambda_{m}-\lambda_{n})t)}{2},$
(S10)

given by the Jensen’s inequality and the right hand side of which is a
function that decreases monotonically with the average eigenvalue gap.
Therefore, the eigenstate witness is also a great indicator to assess the
converged performance of the excited-state VQE (more details can be found in
the Supplemental Material of [74]).

In experiments, the eigenstate witness $r$ can be measured using randomized
measurements by measuring quantum states in $m$ randomly chosen single qubit
basis. We prepare the quantum state $\rho_{\mathrm{reduced}}$ and rotate it
into a random basis with the unitary $V^{(n)}$. Then we measure $s$ samples of
the rotated state $V^{(n)}\rho_{\mathrm{reduced}}V^{(n)^{{\dagger}}}$ in the
computational basis and estimate the probability $P^{(n)}(v_{k})$ of measuring
the computational basis state $v_{k}$ for the basis $n$. This procedure is
repeated for $m$ different measurement bases, and $r$ is then calculated as:

$\displaystyle
r=\sum_{v_{k},v_{q}}(-2)^{-D(v_{k},v_{q})}\sum_{n=1}^{m}P^{(n)}(v_{k})P^{(n)}(v_{q}),$
(S11)

where $D(k,q)$ is the Hamming distance that counts the number of bits that
differ between the computational states $v_{k}$ and $v_{q}$. Because in our
study, we only need to measure the one-qubit density matrix:

$\displaystyle\rho_{\mathrm{reduced}}=\left(\begin{matrix}a&b\\\
b^{*}&1-a\end{matrix}\right),$ (S12)

we can determine the density matrix by the expectation value $\langle
x\rangle$, $\langle y\rangle$ and $\langle z\rangle$ after a simple algebraic
calculation:

$\displaystyle a$ $\displaystyle=$ $\displaystyle\frac{\langle
z\rangle+1}{2},$ (S13) $\displaystyle b$ $\displaystyle=$
$\displaystyle\frac{\langle x\rangle}{2}+\frac{\langle y\rangle}{2}i.$ (S14)

###  .6 F. Quantum noise and Trotter decomposition

In the eigenstate witness measurement case, we have to discuss in terms of
density matrix language when noise is considered. After time evolution
controlled by ancilla qubit, the reduced density matrix can be written as:

$\displaystyle\rho_{\text{reduced}}=\begin{bmatrix}1/2\mathrm{Tr}(\rho)&1/2\mathrm{Tr}(\rho
U^{\dagger})\\\ 1/2\mathrm{Tr}(U\rho)&1/2\mathrm{Tr}(U\rho
U^{\dagger})\end{bmatrix},$ (S15)

where $U=\mathrm{e}^{iHt}$. Considering the noise added by depolarizing
channel, the reduced density matrix is rescaled as:

$\displaystyle\rho_{\text{reduced}}=(1-p_{\mathrm{eff}})\rho_{\text{reduced}}+p_{\mathrm{eff}}I/2,$
(S16)

where $p_{\mathrm{eff}}=33Np$. The coefficient 33 of $p_{\mathrm{eff}}$ is
determined by the number of two-qubit gates in the controlled time evolution
circuit. Constructing an optimal quantum circuit for a general two-qubit gate
requires at most 3 CNOT gates and 15 elementary one-qubit gates [91]. After
adding the control qubit, it requires at most 3 CCNOT gates and 15 two-qubit
gates. Each CCNOT gate can further decompose as 6 CNOT gates [92]. So the
total number of two-qubit gates is 33. If we also consider the effect of
Trotter decomposition, the unitary time evolution matrix can be written as:

$\displaystyle U(t)=\prod_{i}^{N}\mathrm{e}^{iH_{i,i+1}t},$ (S17)

where
$H_{i,i+1}=\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}+V_{0}\sigma^{z}_{i}\sigma^{z}_{i+1}+\frac{h_{i}}{2}\sigma^{z}_{i}+\frac{h_{i+1}}{2}\sigma^{z}_{i+1}$,
$h_{i}=W\cos(2\pi\eta i+\phi)$. In this part, we can also use a PQC with fewer
two-qubit gates to approximate the controlled time evolution and it may
further reduce the effects of noise.

###  .7 G. Details for the real hardware experiments

To extend our method to the real hardware, we need to optimize the parameters
of PQC for the excited-state VQE and the controlled time evolution VQC (see
Fig. S11). The former parameters are optimized as discussed above; the latter
are optimized via the normal VQE classically with following cost function:

$\displaystyle C(\theta)=-\frac{\mathrm{Tr}(U(\theta)V^{\dagger})}{2^{N}},$
(S18)

where $U(\theta)$ is the unitary ansatz of VQC and $N$ is the number of
qubits, $V$ is the controlled time evolution:

$\displaystyle V=\left(\begin{matrix}I&0\\\
0&\mathrm{e}^{iHt}\end{matrix}\right).$ (S19)

Because $U(\theta)V^{\dagger}$ is also a unitary matrix of dimension $2^{N}$,
the cost function has a minimum of $-1$ when $U(\theta)=V$. We can use
$U(\theta)$ as the ansatz to approximate the controlled time evolution in the
real hardware. In our study, we set $\mathrm{D}_{\mathrm{VQE}}=1$ because of
the small system size and $\mathrm{D}_{\mathrm{VQC}}=6$ for a better
approximation (fidelity around 0.97). Only $\mathrm{D}_{\mathrm{VQC}}*N$
($6*4$) two-qubit gates are required for the controlled time evolution VQC
while $33*N$ ($33*4$) two-qubit gates are required for the controlled time
evolution with one time slice Trotter decomposition. Actually, there is a
tradeoff between the depth $\mathrm{D}_{\mathrm{VQC}}$ and the fidelity when
using VQC to approximate the controlled time evolution: deeper depth
$\mathrm{D}_{\mathrm{VQC}}$ can achieve higher fidelity, in the meanwhile,
more resources are required. As suggested by the results obtained from noisy
simulation and the theoretical understanding of eigenstate witness, we don’t
need perfect controlled time evolution to discriminate the two phases. The VQC
with relatively high fidelity is enough to distinguish the MBL and the thermal
phases. Therefore, for larger systems, we can still utilize a VQC to
approximate the controlled time evolution with much fewer two-qubit gates than
those required in Trotter decomposition.

We set $W=8.0$ for the MBL phase and $W=1.5$ for the thermal phase, and we
choose the evolved time $t=0.15/W$. The parameters of the excited-state VQE
and VQC are both determined by the best of $20$ independent VQE trials. In
Fig. 6, the results of the noisy simulator are the average of $100*8192$
independent measurement shots and the results of real hardware are the average
of $50*8192$ independent measurement shots.

Figure S11: VQC using hardware-efficient ansatz to approximate controlled time
evolution. The two-qubit entangling gate is $\mathrm{CZ}$ gate.