# Efficient criteria of quantumness for a large system of qubits

Shohei Watabe1,2, Michael Zach Serikow3,4, Shiro Kawabata5, and Alexandre
Zagoskin6 1 Department of Physics, Faculty of Science Division I, Tokyo
University of Science, Shinjuku, Tokyo 162-8601, Japan  2 Division of Nano-
quantum Information Science and Technology, Research Institute for Science and
Technology, Tokyo University of Science, Shinjuku, Tokyo 162-8601, Japan  3
Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA  4
Notre Dame Institute for Advanced Study, University of Notre Dame, IN 46556,
USA  5 Research Center for Emerging Computing Technologies (RCECT), National
Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono,
Tsukuba, Ibaraki 305-8568, Japan.  6 Department of Physics, Loughborough
University, Loughborough LE11 3TU, UK

###### Abstract

In order to model and evaluate large-scale quantum systems, e.g. quantum
computer and quantum annealer, it is necessary to quantify the “quantumness”
of such systems. In this paper, we discuss the dimensionless combinations of
basic parameters of large, partially quantum coherent systems, which could be
used to characterize their degree of quantumness. Based on model analytical
and numerical calculations, we suggest one such number for a system of qubits
undergoing adiabatic evolution. Applying it to the case of D-Wave One
superconducting quantum annealing device, we find that its operation as
described falls well within the quantum domain.

##  I Introduction

One of the key obstacles in the way to the full development of quantum
technologies 2.0 is the same circumstance which stimulated their development
in the first place: the fundamental impossibility of an efficient simulation
of a large enough, quantum coherent structure with classical means. In
practice “large enough” turned out to be a system comprising a hundred or so
quantum bits, which is still too small to form a quantum computer capable of
simulating other “large enough” quantum systems. On the other hand, artificial
quantum coherent systems comprising thousands of qubits [1] are being
fabricated and even successfully used, like commercial quantum annealers [2,
3]. This “quantum capacity gap” [4] needs to be bridged, in order to allow a
systematic progress towards the development of the full potential of quantum
technologies 2.0, such as universal quantum computers.

The impossibility of an efficient classical simulation of a large quantum
system is not absolute, in the sense that it concerns the simulation of an
arbitrary evolution of such a system, whereby its state vector can reach all
of its (exponentially high-dimensional) Hilbert space and has potentially
infinite time to do so. The Margolus-Levitin theorem and its generalizations
[5, 6, 7, 8, 9, 10] put a limit on the speed of such evolution, thus
restricting the accessible part of the Hilbert space for any finite time
interval. This agrees with a proof [11] that the manifold of all quantum many-
body states that can be generated by arbitrary time-dependent local
Hamiltonians in a time that scales polynomially in the system size occupies an
exponentially small volume in its Hilbert space. (This is a literally correct
statement, since the Hilbert space of a system of qubits is a finite-
dimensional complex projective space; that is, it is compact and, moreover, it
has a unitary invariant Fubini-Study metric [12].) Numerical and analytical
studies also indicate that the number of independent constraints describing
quantum evolution may be much less than the dimensionality of the Hilbert
space [13]. It is therefore reasonable to suggest that a “general case”
evolution of a large quantum coherent quantum structure can be characterized
by a non-exponentially large set of dimensionless parameters, which correspond
to qualitatively different regimes of evolution of this structure. Our recent
numerical simulations indicate the existence of such regimes in a set of
qubits with pumping and dissipation [14].

Such dimensionless parameters, if exist, will be combinations of fundamental
physical constants and parameters, which characterize the system. (We will
only need the Planck constant, since the speed of light and gravity constant
are not relevant for the currently feasible devices). For example, a system of
qubits is generally described by a quantum Ising Hamiltonian,

$\displaystyle\hat{H}(t)=\sum_{i\neq
j}J_{ij}(t)\sigma_{i}^{z}\sigma_{j}^{z}-\frac{1}{2}\sum_{j}\left(h_{j}(t)\sigma_{j}^{z}+\Delta_{j}(t)\sigma_{j}^{x}\right),$
(1)

and a set of Lindblad operators responsible for dephasing and relaxation of
separate qubits, with characteristic times, respectively, $t_{\phi}$ and
$t_{r}$, and their combination, the decoherence time $t_{D}$. In the case of
an adiabatic quantum processor [2, 3], the time dependence of the Hamiltonian
parameters $J_{ij},h_{j}$ and $\Delta_{j}$ (except that induced by the ambient
noise) is determined by that of the adiabatic parameter $\lambda(t)$ (in the
simplest case, by the time of adiabatic evolution, $t_{f}$). Then, the
dimensionless characteristics of the system should be the combinations of
$\hbar$ with the following quantities:

  1. 1.

Dimensionless: $N$ (number of qubits); $\langle z\rangle$ (connectivity of the
network: average number of couplings per qubit); $\langle\delta z^{2}\rangle$
(its dispersion); $\langle z_{i}z_{j}\rangle$ (its correlation function);
$\dots$;

  2. 2.

Powers of energy: $\langle E\rangle$ (average qubit excitation energy);
$\langle\delta E^{2}\rangle$; $\langle E_{i}E_{j}\rangle$; $\langle J\rangle$
(average coupling strength); $\langle\delta J^{2}\rangle$; $\langle
J_{i}J_{k}\rangle$; $k_{B}T$; $\dots$;

  3. 3.

Powers of time: $t_{\phi}$; $t_{r}$; $t_{D}$; $t_{f}$; $\dot{\lambda}\sim
1/t_{f}$ (speed of evolution); $\ddot{\lambda}$; $\dots$;

and such additional parameters as, e.g., the spectral density of ambient noise
$S_{A}(f)$. Note that all these parameters can be efficiently obtained by
either direct measurements or straightforward calculations.

The field of several dozens of independent dimensionless combinations of the
above parameters is narrowed for a particular quantum system and the mode of
its operation. Here we concentrate on adiabatic quantum computing (quantum
annealing). The Hamiltonian of a quantum system is here manipulated in such a
way that its ground state changes from the easily accessible one to the one
encoding a solution to the desired problem (and presumably having a very
complex structure, so that reaching it by annealing would be improbable). In
the case of a slow enough evolution of the Hamiltonian $H(t)$, the system
initialized in the ground state of $H(0)$ will evolve into the ground state of
$H(t_{f})$ by virtue of the adiabatic theorem [15]. If the system is totally
insulated, the quantum speed and accessibility theorems [5, 6, 7, 8, 9, 10,
11] do not put fundamental constraints on an adiabatic quantum computer.
Nevertheless, in a realistic case, the computation time is limited by
interactions with the outside world leading to nonunitarity [16], and the
question arises whether the evolution from the initial to the desired final
state of the system is possible. As we see, the question of accessibility of
different regions of the Hilbert space is especially relevant in this case.

An intriguing twist is added by the fact that the operation of D-Wave
processors demonstrated what looked convincingly like quantum annealing [17,
18] despite the large discrepancy between the adiabatic evolution time,
$t_{f}$ (microseconds [17]), and the qubit decoherence time, $t_{D}\ll t_{f}$
(tens of nanoseconds, [19]), in the absence of any quantum error correction.
On the second thought, it is not so surprising. The quantum state of an
adiabatic quantum computer evolves starting at $t=0$ from a factorized state,
$|{\rm in}\rangle$. The computation is successful, if at the time $t=t_{f}$
there is a sufficient ratio of quantum trajectories ending in the state $|{\rm
out}\rangle$, which is also factorized by design. Decoherence tends to disrupt
quantum correlations between different qubits, and thus constrain the
trajectories to partially factorized submanifolds of the Hilbert space.
Nevertheless the success does not necessarily require that these trajectories
pass through globally entangled states, and thus certain degree of decoherence
may not necessarily make the proper operation of the device impossible.

While it cannot be predicted whether the evolution of a given quantum system
can take it from the given initial to the desired final state, the average of
the maximal distance between _some_ initial and _some_ final state of the
system, for given values of $t_{f}$ and other system parameters, may serve as
a heuristic indicator of success. This distance can be naturally determined
via the Fubini-Study metric [12], in which the distance $s(\phi,\psi)$ between
states $|\phi\rangle$ and $|\psi\rangle$ is given by

$\displaystyle\cos
s(\phi,\psi)=\sqrt{\frac{\langle\phi|\psi\rangle\langle\psi|\phi\rangle}{\langle\phi|\phi\rangle\langle\psi|\psi\rangle}}.$
(2)

The maximal Fubini-Study distance in the Hilbert space is $\pi/2$, the
distance between mutually orthogonal states. We will therefore choose the
quantity

$\displaystyle{\cal S}=\frac{2}{\pi}\sqrt{\overline{s^{2}(\phi,\psi)}}$ (3)

as an _ad hoc_ parameter, which characterizes the ability of an adiabatically
evolving quantum device to reach its desired quantum state. The bar denotes
the averaging over all initial states and over all quantum trajectories
accessible to the system, which connect them to the final states maximally
removed from them (in terms of the Fubini-Study distance).

##  II Random walk model: Heuristic treatment

The evolution of the state vector of a quantum system can be modeled by a
series of random collapses to one of its instantaneous eigenstates at the
moments $t_{1},t_{2},\dots$, and unitary evolutions under the Hamiltonian
$\hat{H}(t)$ between these moments. It was very tempting to relate the finite
probability of success of an adiabatic quantum computer to the percolation
transition in its Hilbert space when covered by randomly placed hyperspheres
of the radius, determined by the interval $\Delta t$ between the collapses. A
formation of the equivalent to the infinite cluster [20] would guarantee the
finite probability of accessing the desired final state. Unfortunately, we
have concluded that the approach does not scale, though the results by
themselves are instructive (see Section IV).

The next simplest model is a random walk in the Hilbert space of the system,
with the step (Fubini-Study) lengths $\Delta s_{j}$ dependent on $\Delta t$.
For small time intervals of unitary evolution $|\phi\rangle=e^{-i\hat{H}\Delta
t}|\psi\rangle$, the distance $s(\phi,\psi)$ is given by (see Appendix A)

$\displaystyle s\approx\frac{\Delta t}{\hbar}\sigma_{\hat{H}}(\psi).$ (4)

Here,
$\sigma_{\hat{H}}^{2}(\psi)\equiv\langle\hat{H}^{2}\rangle_{\psi}-\langle\hat{H}\rangle_{\psi}^{2}$
is the energy dispersion during the unitary evolution, where
$\langle\,\cdot\,\rangle_{\psi}$ is the expectation value for a quantum state
$|\psi\rangle$, i.e.,
$\langle\,\cdot\,\rangle_{\psi}\equiv\langle\psi|\cdot|\psi\rangle$. This
expression is intuitively plausible: energy and time are the only parameters
in the problem, and in case of zero energy variance the system would be in an
eigenstate of the Hamiltonian and remain in it, barring degeneracies.

Substituting here the maximal possible value of $s=\pi/2$ (which is, strictly
speaking, not allowed), we find

$\displaystyle\Delta t_{\rm
max}=\frac{\pi}{2}\frac{\hbar}{\sigma_{\hat{H}}(\psi)}.$ (5)

This reproduces the rigorous Mandelstam-Tamm expression for the minimal time
necessary to evolve from an initial state to a state orthogonal to it [5].

Now consider the random walk of $M\gg 1$ steps of identical duration $\Delta
t$, controlled by independently distributed random Hamiltonians $\hat{H}_{j}$
($j=1,2,\dots M$) (but still assuming that the total displacement remains
small). Using the same approximation as before, the span of this random walk
(i.e., the Fubini-Study distance between its initial point $|\psi\rangle$ and
final point $|\phi\rangle=\prod_{j}\exp{(-i\Delta
t\hat{H}_{j}/\hbar)}|\psi\rangle$) is thus found directly (see Appendix B):

$\displaystyle s(\phi,\psi)\approx\frac{\Delta
t}{\hbar}\sigma_{\sum_{j}\hat{H}_{j}}(\psi),$ (6)

where

$\displaystyle\sigma_{\sum_{j}\hat{H}_{j}}(\psi)=$
$\displaystyle\biggl{\langle}\sum_{j}\hat{H}_{j}^{2}\biggr{\rangle}_{\psi}-\biggl{\langle}\sum_{j}\hat{H}_{j}\biggr{\rangle}_{\psi}^{2}$
(7) $\displaystyle=$
$\displaystyle\sum_{\chi\neq\psi}\sum_{jk}\langle\psi|\hat{H}_{j}|\chi\rangle\langle\chi|\hat{H}_{k}|\psi\rangle.$
(8)

Making a further simplification, assume that all random steps have the same
length $\Delta S=\overline{\sigma_{\hat{H}}}\Delta t/\hbar$. Then the r.m.s.
of the span is

$\displaystyle\sqrt{\overline{s(\phi,\psi)^{2}}}\approx\Delta
S\frac{\sqrt{\overline{\sigma_{\sum_{j}\hat{H}_{j}}^{2}(\psi)}}}{\overline{\sigma_{\hat{H}}(\psi)}},$
(9)

where the bar average is over random choice of $H_{j}$-Hamiltonians and
initial states $|\psi\rangle$.

After averaging over random Hamiltonians $\hat{H}_{j}$ and states
$|\psi\rangle$, which gives
$\overline{\sum_{j}\langle\psi|\hat{H}_{j}|\chi\rangle}\approx 0$, the only
term surviving will be
$\overline{\sigma_{\sum\hat{H}_{j}}^{2}(\psi)}=\sum_{\chi\neq\psi}\sum_{j}\overline{|\langle\psi|\hat{H}_{j}|\chi\rangle|^{2}}$.
Since the averaged value with respect to random $\hat{H}_{j}$ is
$j$-independent, we have

$\displaystyle\overline{\sigma_{\sum\hat{H}_{j}}^{2}(\psi)}=$
$\displaystyle\sum_{j=1}^{M}\sum_{\chi\neq\psi}\overline{|\langle\psi|\hat{H}|\chi\rangle|^{2}}$
(10) $\displaystyle\approx$ $\displaystyle
Mf(N)\overline{\sigma_{\hat{H}}(\psi)}^{2},$ (11)

where $f(N)$ is some function of the dimension of the Hilbert space $D=2^{N}$.
Then we obtain from Eq. (3)

$\displaystyle{\cal S}\approx\frac{2}{\pi}\Delta SM^{1/2}f(N)^{1/2}.$ (12)

In particular, the condition ${\cal S}=1$ yields the relation between $M$ and
$\Delta S$, for which the random walk is likely to connect mutually orthogonal
states in the Hilbert space (a _maximal random walk_), and thus the quantum
adiabatic operation of the system we model should be possible. This “critical
value” of $\Delta S$ for a given $M,N$ is

$\displaystyle\Delta S=\frac{\pi}{2}M^{-1/2}f(N)^{-1/2}.$ (13)

Figure 1:  Random walk in the Hilbert space. (a) Schematics of states walking
randomly in the Hilbert space. The points represent states. The state walks
randomly from $|\Psi_{t=0}\rangle$ to $|\Psi_{t=M}\rangle$, where the stride
of each step is given by the Fubini-Study distance $\Delta S$. (b) Schematic
of the random walk in the Hilbert space in the Bloch sphere. The red points
indicate the states, which walks with the random direction with a fixed stride
$\Delta S$. The random walker on the Bloch sphere starts from the initial
state $|0\rangle$ at the north pole, and finally reaches the state $|1\rangle$
at the south pole. (c) For considering the uniform random walk in the Hilbert
space, we generate a random vector in the local orthogonal coordinate (spanned
by the bases ${\bf e}_{1,2}^{\prime}$, where ${\bf e}_{1}^{\prime}\cdot{\bf
e}_{2}^{\prime}=0$, for example) from the curved parameter space
${\bm{\theta}}$ (spanned by the bases ${\bf e}_{1,2}$, where ${\bf
e}_{1}\cdot{\bf e}_{2}\neq 0$). (d) The critical Fubini-Study distance $\Delta
S$ for the random walk in the in the Hilbert space, as a function of the
number of random steps $M$. We consider $N$-qubits states. The data points are
the averaged value for $500$ trials. Solid lines are fitted lines with $\Delta
S=(\pi/2)AM^{-B}$. The fitting is performed for $M\geq 3$. (e) The fitting
factor $A$ and (f) the fitting factor $B$ as a function of the number of
qubits $N$. The solid lines represent fitting functions, where we used
$A=\alpha_{A}N^{\beta_{A}}$ and $B=\alpha_{B}N^{\beta_{B}}$ for $N=1$–$6$.

##  III Random walk model: Numerical approach

The goal of our numerical simulations is to determine the relation between the
Fubini-Study length of single step, $\Delta S$, and the number of steps, $M$,
of a maximal random walk in the $D=2^{N}$-dimensional Hilbert space (see Figs.
1 (a) and (b)).

First, we randomly prepare an initial state $|\Psi_{t=0}\rangle$ in the
$D$-dimensional Hilbert space. The state $|\Psi_{t}({\bm{\theta}})\rangle$ on
step $t$ is parameterized by the $2(D-1)$-dimensional real vector

$\displaystyle{\bm{\theta}}=(\theta_{1},\theta_{2},\cdots,\theta_{2(D-1)})^{\rm
T},$ (14)

so that $|\Psi_{t+1}({\bm{\theta}}+d{\bm{\theta}})\rangle$.

The Fubini-Study step length is then

$\displaystyle\Delta
S=s(\Psi_{t}({\bm{\theta}}),\Psi_{t+1}({\bm{\theta}}+d{\bm{\theta}}))=\sum_{ij}g_{ij}d\theta_{i}d\theta_{j},$
(15)

where the Fubini-Study metric
$g_{ij}=\langle\partial_{\theta_{i}}\Psi({\bm{\theta}})|\partial_{\theta_{j}}\Psi({\bm{\theta}})\rangle-\langle\partial_{\theta_{i}}\Psi({\bm{\theta}})|\Psi({\bm{\theta}})\rangle\langle\Psi({\bm{\theta}})|\partial_{\theta_{j}}\Psi({\bm{\theta}})\rangle$,
and
$|\partial_{\theta_{i}}\Psi({\bm{\theta}})\rangle\equiv\partial|\Psi({\bm{\theta}})\rangle/\partial\theta_{i}$.
We generate the parameter $d{\bm{\theta}}$ such that the direction is
uniformly random in the local orthogonal space (see Figs. 1 (c)), and the
Fubini-Study distance $\sum_{ij}g_{ij}d\theta_{i}d\theta_{j}$ satisfies a
given fixed value $\Delta S$. (See Appendix C.) We will keep the state vectors
normalized. The span of the random walk is then

$\displaystyle s(\Psi_{0},\Psi_{t})=\arccos(\langle\Psi_{0}|\Psi_{t}\rangle).$
(16)

In our simulations we fix the number of steps $M$ and search for such a value
of $\Delta S$ (“critial value”), that the final state $|\Psi_{t=M}\rangle$
satisfies the condition $|\pi/2-s(\Psi_{0},\Psi_{M})|\leq\epsilon$ for some
small $\epsilon$.

Figure 1 (d) shows the critical value of $\Delta S$ as a function of $M$.
Fitting yields

$\displaystyle\Delta S=\frac{\pi}{2}A(N)M^{-B(N)},$ (17)

where $A(N)=0.309(8)N^{0.76(2)}$ and $B(N)=0.4727(8)N^{0.020(1)}$ (see Figs 1
(e) and (f)). In these numerical simulations, we used
$\epsilon=\pi/2-\cos^{-1}(1/\sqrt{D})$, where $\epsilon\approx 1/\sqrt{D}$ in
the huge Hilbert space. We see that the critical value of $\Delta S$ scales
almost as a power of $M$, since $B(N)$ shows a very weak power law dependence
on the number of qubits.

Comparing Eqs.(17) and (13), we see that the dependence on $M$ in our
heuristic and numerical approaches is almost the same, while for the function
$f(N)$ we find $f(N)\approx 10.5\>N^{-3/2}$.

##  IV Accessibility index, ”quantumness” criterion and comparison to
experiment

The quantities $\Delta S$ and $M$ do not have a direct experimental
significance. In our approximate treatment we can relate them instead to the
decoherence time (during which the unitary evolution takes place), $t_{D}$,
and the total time of adiabatic evolution, $t_{f}$, via

$\displaystyle M=\frac{t_{f}}{t_{D}};\>\>\Delta S\leq\frac{t_{D}}{\hbar}J,$
(18)

where $J$ is the typical coupling between qubits, so that $J$ is a reasonable
measure of the uncertainty of the $N$-qubit system’s energy during the
adiabatic evolution (here we use the Mandelstam-Tamm expression). Then from
(17) we find, that the necessary “quantumness” criterion for the adiabatic
evolution is

$\displaystyle t_{D}\frac{2\pi
J}{h}>\frac{\pi}{2}CN^{3/4}\left(\frac{t_{D}}{t_{f}}\right)^{1/2},$ (19)

where $C\approx 0.31$. Given all the approximations we have made, we can as
well take $C=1$. Then the ”quantumness” condition can be written as

${\cal A}>1,$ (20)

where the _accessibility index_ for a system of $N$ qubits with average
coupling strength $J$ is

${\cal A}\equiv 4\frac{J}{h}\frac{\sqrt{t_{f}t_{D}}}{N^{3/4}}.$ (21)

Applying this criterion to the operation of the D-Wave processor described in
Refs. [19, 17], with $t_{D}\sim 10$ ns, $t_{f}\sim 5\mu$s, $N\sim 100$ and
$J/h\sim 5$GHz, we see that

$\displaystyle{\cal A}\approx 10^{2}\gg 1,$ (22)

and the necessary “quantumness” condition was satisfied. This indicates that
the results of Refs. [19, 17] are consistent with quantum annealing. From Eq.
(20), we can evaluate the maximal size of a quantum processor for which the
”quantumness” condition holds, other things being the same as in [19, 17]:
$N_{\rm max}\equiv[4J/h]^{4/3}(t_{f}t_{D})^{2/3}\approx 7\times 10^{4}$.

Note that the condition ${\cal A}>1$ does not guarantee the successful
operation of a quantum processor, i.e., it reaching the desired final state.
It is rather the necessary, but not sufficient, condition of success.

##  V Discussion: Percolation model

Figure 2:  Percolation in the Hilbert space. (a) Schematics of connection
between states. The points represent states generated randomly, and the
diameter of the circle shows the distance $\Delta S$. Red/blue shows
connected/disconnected states. (b) Schematic of the percolation on the Bloch
sphere. The red/blue points indicate the random states generated uniformly on
the Bloch sphere, where the red/blue points can/cannot be organized into the
cluster maximally spans on the Bloch sphere. The diameter of the red and blue
circles is $\Delta S$. The red cluster spans from the state $|0\rangle$ to the
orthogonal state $|1\rangle$. (c) The Fubini-Study threshold distance $\Delta
S$ for connecting two states as a function of the number of isotropic random
states $M$ in the Hilbert space for $N$-qubits, where the cluster can expand
the maximum Fubini-Study distance $\pi/2$. The data points are the averaged
value for $100$ samples, and solid lines are fitted lines with $\Delta
S=(\pi/2)AM^{-B}$ for $M=2$–$200$. (d) Schematics of the Fubini-Study distance
between two isotropically random states in the high-dimensional Hilbert space.
These two states are almost orthogonal. (e) The fitting factor $A$ and (f) the
fitting factor $B$ as a function of the number of qubits $N$ and the dimension
of the Hilbert space: $D=2^{N}$. The solid lines represent fitting functions,
where we used $A=\gamma_{A}-\alpha_{A}D^{-\beta_{A}}$ and
$B=\alpha_{B}D^{-\beta_{B}}$ for $N=7$–$14$. The inset in (e) is a log-log
plot of $\gamma_{A}-A$ as a function of $D$.

Here, for completeness, we consider a percolation model of adiabatic state
evolution in a partially coherent system. The simplest case would be a
continuous percolation through randomly placed hyperspheres. Due to
compactness, and actually finite diameter of the Hilbert space [12], the
standard definition of a percolation transition (appearance of the infinite
cluster [20]) does not apply. Instead we introduced the notion of “maximal
span cluster”, MSC, that is, a cluster which contains at least two points
separated by the maximal Bures angle, $\pi$.

However, the percolation transition is ruled out as a useful model of
evolution of large, partially coherent quantum structures, because two
uniformly generated random states $|\Psi\rangle$ and $|\Psi^{\prime}\rangle$
become almost orthogonal as the dimensionality of the Hilbert space grows. The
concentration inequality for the absolute square of the inner product between
these two states (see Appendix D) is

$\displaystyle
P\left(\left|\langle\Psi|\Psi^{\prime}\rangle|^{2}-\frac{1}{D}\right|\geq\epsilon\right)\leq
4\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right],$
(23)

where $|\langle\Psi|\Psi^{\prime}\rangle|^{2}$ is almost the inverse of the
dimension $D=2^{N}$, i.e. exponentially small for the number of qubits
$|\langle\Psi|\Psi^{\prime}\rangle|^{2}={\mathcal{O}}(2^{-N})$. Taking
$\epsilon=1/D$, we have $P_{\rm L.B.}\leq
P\left(0<|\langle\Psi|\Psi^{\prime}\rangle|^{2}<2/2^{N}\right)\leq 1$, where
$P_{\rm L.B.}\equiv 1-4\exp\left(-2^{N}/16\right)$. For example, the lower
bound $P_{\rm L.B.}$ is $0.92\cdots$ for $6$-qubits, $0.998\cdots$ for
$7$-qubits, and $0.9999995\cdots$ for $8$-qubits. The probability where the
absolute square of this inner product ranges from $0$ to $2/2^{N}$ is double-
exponentially approaching to unity for large $N$.

We formulated and investigated numerically the continuum percolation
transition in a finite-dimensional Hilbert space randomly covered by
hyperspheres (see Appendix E for details). We first randomly prepare
$M$-states $|\Psi_{n}\rangle$ for $n=1,2,\cdots M$, which are uniformly
distributed in the Hilbert space. We then create clusters by connecting two
states $|\Psi_{m,n}\rangle$ with the condition where the Fubini–Study distance
$s(\Psi_{m},\Psi_{n})=\cos^{-1}(|\langle\Psi_{m}|\Psi_{n}\rangle|)$ is shorter
than a certain threshold $\Delta S$, i.e., $s(\Psi_{m},\Psi_{n})\leq\Delta S$
(see Fig. 2 (a)). Let $\alpha$ be the label of the cluster. By comparing the
distances between all the states in the $\alpha$-th cluster, we can find the
maximum distance $L_{\alpha}={\rm max}(s(\Psi_{m},\Psi_{n}))$, where
$\Psi_{m,n}$ are the states that belong to the $\alpha$-th cluster. If we find
$|\pi/2-L_{\alpha}|\leq\epsilon$ for some small $\epsilon$, we regard that the
$\alpha$-th cluster spans maximally in the Hilbert space (see Fig. 2 (b)).

Figure 2 (c) plots the critical Fubini-Study distance $\Delta S$ as a function
of the number of the random states $M$ generated uniformly in the Hilbert
space. Here, the critical Fubini-Study distance is determined by the threshold
of the distance $\Delta S$, where a cluster starts to span maximally with the
condition $|\pi/2-L_{\alpha}|\leq\epsilon$. In this numerical simulation, we
take $\epsilon=\Delta S$. The scaling of the critical Fubini-Study distance
$\Delta S$ can be given by $\Delta S=(\pi/2)AM^{-B}$, where
$A(D)=1.0005(5)-0.33(2)D^{-0.47(2)}$ and $B(D)=0.182(4)D^{-0.522(4)}$, with
the dimension of the Hilbert space $D=2^{N}$. A numerical simulation shows
that the exponent scales to zero as a power of the dimension of the Hilbert
space, which unfortunately rules out the percolation transition as a candidate
for the adiabatic accessibility index. In a system of many qubits, we need
almost the maximum value of the distance $\Delta S\simeq\pi/2$ to make the
cluster span maximal irrespective of the number of random states $M$, which is
consistent with the insight of the concentration inequality in Eq. (23).

##  VI Conclusions

We have investigated the generic behavior of a partially coherent system of
qubits undergoing adiabatic evolution. Our aim was finding a convenient
dimensionless parameter, which could characterize the degree of “quantumness”
of our system.

One such promising parameter related to the percolation approach in Hilbert
space was found to be not scalable.

Basing our analysis instead on a random walk model of quantum evolution in the
Hilbert space, we found a parameter allowing to evaluate the likelihood of a
successful quantum transition between the initial and desired final states of
the system. It is the _accessibility index_ , expressed through the qubit
decoherence time, time of evolution and the qubit coupling strength. Applying
it to the case of 128-qubit D-Wave processors, we found that their evolution
was consistent with quantum adiabatic transitions despite the qubit
decoherence time being much smaller than the evolution time.

###### Acknowledgements.

This paper is partly based on results obtained from a project, JPNP16007,
commissioned by the New Energy and Industrial Technology Development
Organization (NEDO), Japan. S.W. was supported by Nanotech CUPAL, National
Institute of Advanced Industrial Science and Technology (AIST). A.Z. was
supported by NDIAS Residential Fellowship.

## References

  * Kakuyanagi  _et al._ [2016] K. Kakuyanagi, Y. Matsuzaki, C. Déprez, H. Toida, K. Semba, H. Yamaguchi, W. J. Munro, and S. Saito, Observation of collective coupling between an engineered ensemble of macroscopic artificial atoms and a superconducting resonator, Phys. Rev. Lett. 117, 210503 (2016).
  * [2] Introduction to the D-Wave Quantum Hardware, https://www.dwavesys.com/tutorials/background-reading-series/introduction-d-wave-quantum-hardware.
  * Bunyk  _et al._ [2014] P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, and J. Whittaker, Architectural Considerations in the Design of a Superconducting Quantum Annealing Processor, IEEE Transactions on Applied Superconductivity 24, 1 (2014).
  * Zagoskin [2014] A. Zagoskin, The Grand Challenge of Quantum Computing: Bridging the Capacity Gap, Frontiers in ICT 1, 2 (2014).
  * Mandelstam and Tamm [1945] L. Mandelstam and I. Tamm, The uncertainty relation between energy and time in nonrelativistic quantum mechanics, Izv. Akad. Nauk SSSR ser. fiz. 9, 122 (1945).
  * Margolus and Levitin [1998] N. Margolus and L. B. Levitin, The maximum speed of dynamical evolution, Physica D: Nonlinear Phenomena 120, 188 (1998).
  * Zwierz  _et al._ [2012] M. Zwierz, C. A. Pérez-Delgado, and P. Kok, Ultimate limits to quantum metrology and the meaning of the Heisenberg limit, Physical Review A 85, 1 (2012).
  * Deffner and Lutz [2013a] S. Deffner and E. Lutz, Energy–time uncertainty relation for driven quantum systems, Journal of Physics A: Mathematical and Theoretical 46, 335302 (2013a).
  * Deffner and Lutz [2013b] S. Deffner and E. Lutz, Quantum Speed Limit for Non-Markovian Dynamics, Physical Review Letters 111, 249 (2013b).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: from Heisenberg’s uncertainty principle to optimal quantum control, Journal of Physics A: Mathematical and Theoretical 50, 453001 (2017).
  * Poulin  _et al._ [2011] D. Poulin, A. Qarry, R. Somma, and F. Verstraete, Quantum simulation of time-dependent hamiltonians and the convenient illusion of hilbert space, Phys. Rev. Lett. 106, 170501 (2011).
  * Brody and Hughston [2001] D. C. Brody and L. P. Hughston, Geometric quantum mechanics, Journal of Geometry and Physics 38, 19 (2001).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Andreev  _et al._ [2021] A. V. Andreev, A. G. Balanov, T. M. Fromhold, M. T. Greenaway, A. E. Hramov, W. Li, V. V. Makarov, and A. M. Zagoskin, Emergence and control of complex behaviors in driven systems of interacting qubits with dissipation, npj Quantum Information 7, 1 (2021).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
  * Sarandy and Lidar [2005] M. S. Sarandy and D. A. Lidar, Adiabatic quantum computation in open systems, Phys. Rev. Lett. 95, 250503 (2005).
  * Boixo  _et al._ [2014] S. Boixo, T. F. Ronnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nature Physics 10, 218 (2014).
  * Albash  _et al._ [2015] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, Consistency tests of classical and quantum models for a quantum annealer, Phys. Rev. A 91, 042314 (2015).
  * Boixo  _et al._ [2013] S. Boixo, T. Albash, F. Spedalieri, N. Chancellor, and D. Lidar., Experimental signature of programmable quantum annealing, Nat. Comm. 4, 3067 (2013).
  * Stauffer and Aharony [1992] D. Stauffer and A. Aharony,  _Introduction To Percolation Theory: Second Edition_ ((London: Taylor & Francis), 1992).
  * Hall  _et al._ [2005] P. Hall, J. S. Marron, and A. Neeman, Geometric representation of high dimension, low sample size data, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67, 427 (2005).
  * Wainwright [2019] M. J. Wainwright,  _High-Dimensional Statistics: A Non-Asymptotic Viewpoint_, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 2019).

##  Appendix A Derivation of Eq. (4)

In this appendix, we derive Eq. (4). From Eq. (2), the cosine distance between
$|\phi\rangle$ and $|\psi\rangle$ is given by $\cos
s(\phi,\psi)=\sqrt{\langle\phi|\psi\rangle\langle\psi|\phi\rangle/(\langle\phi|\phi\rangle\langle\psi|\psi\rangle)}$.
For small time intervals of unitary evolution, we have

$\displaystyle|\phi\rangle=$ $\displaystyle e^{-i\hat{H}\Delta
t/\hbar}|\psi\rangle$ (24) $\displaystyle=$
$\displaystyle\left[1-\frac{i}{\hbar}\hat{H}\Delta
t-\frac{1}{2\hbar^{2}}\hat{H}^{2}(\Delta t)^{2}\right]|\psi\rangle+o((\Delta
t)^{2}).$ (25)

Here,

$\displaystyle\langle\phi|\psi\rangle\langle\psi|\phi\rangle=$
$\displaystyle\left|\langle\psi|\left[1-\frac{i}{\hbar}\hat{H}\Delta
t-\frac{1}{2\hbar^{2}}\hat{H}^{2}(\Delta
t)^{2}\right]|\psi\rangle\right|^{2}+o((\Delta t)^{2})$ (26) $\displaystyle=$
$\displaystyle\left|1-\frac{i\Delta
t}{\hbar}\langle\hat{H}\rangle_{\psi}-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\langle\hat{H}^{2}\rangle_{\psi}\right|^{2}+o((\Delta
t)^{2})$ (27) $\displaystyle=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{\hbar^{2}}\langle\hat{H}^{2}\rangle_{\psi}+\frac{(\Delta
t)^{2}}{\hbar^{2}}\langle\hat{H}\rangle_{\psi}^{2}+o((\Delta t)^{2})$ (28)
$\displaystyle=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{\hbar^{2}}\sigma_{\hat{H}}^{2}(\psi)+o((\Delta t)^{2}).$ (29)

where
$\sigma_{\hat{H}}^{2}(\psi)\equiv\langle\hat{H}^{2}\rangle_{\psi}-\langle\hat{H}\rangle_{\psi}^{2}$.
On the other hand, $\langle\psi|\psi\rangle=1$ and

$\displaystyle\langle\phi|\phi\rangle=$
$\displaystyle\langle\psi|\left[1+\frac{i}{\hbar}\hat{H}\Delta
t-\frac{1}{2\hbar^{2}}\hat{H}^{2}(\Delta
t)^{2}\right]\left[1-\frac{i}{\hbar}\hat{H}\Delta
t-\frac{1}{2\hbar^{2}}\hat{H}^{2}(\Delta t)^{2}\right]|\psi\rangle+o((\Delta
t)^{2})$ (30) $\displaystyle=$
$\displaystyle\langle\psi|\left[1+\frac{1}{\hbar^{2}}\hat{H}^{2}(\Delta
t)^{2}-2\frac{1}{2\hbar^{2}}\hat{H}^{2}(\Delta
t)^{2}\right]|\psi\rangle+o((\Delta t)^{2})$ (31) $\displaystyle=$
$\displaystyle 1+o((\Delta t)^{2}).$ (32)

We can write the right hand side of Eq. (2) as

$\displaystyle\sqrt{\frac{\langle\phi|\psi\rangle\langle\psi|\phi\rangle}{\langle\phi|\phi\rangle\langle\psi|\psi\rangle}}=$
$\displaystyle\sqrt{1-\frac{(\Delta
t)^{2}}{\hbar^{2}}\sigma_{H}^{2}(\psi)}+o((\Delta t)^{2})$ (33)
$\displaystyle=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sigma_{H}^{2}(\psi)+o((\Delta t)^{2}),$ (34)

Since $s=O(\Delta t)$, we have

$\displaystyle\cos s(\phi,\psi)=1-\frac{1}{2}s^{2}+o((\Delta t)^{2}).$ (35)

Therefore using Eqs (2), (34), and (35), we can obtain Eq. (4).

##  Appendix B Derivation of Eq. (6)

Using the same approximation as in Appendix A, we obtain

$\displaystyle|\phi\rangle=$ $\displaystyle\prod_{j}\exp\left(\frac{-i\Delta
t}{\hbar}\hat{H}_{j}\right)|\psi\rangle$ (36) $\displaystyle\approx$
$\displaystyle\prod_{j}\left(1-\frac{i\Delta
t}{\hbar}\hat{H}_{j}-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\hat{H}_{j}^{2}\right)|\psi\rangle$ (37)
$\displaystyle\approx$ $\displaystyle\left(1-\frac{i\Delta
t}{\hbar}\sum_{j}\hat{H}_{j}-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j}\hat{H}_{j}^{2}-\frac{(\Delta
t)^{2}}{\hbar^{2}}\sum_{j>k}\hat{H}_{j}\hat{H}_{k}\right)|\psi\rangle.$ (38)

Here, we have

$\displaystyle|\langle\psi|\phi\rangle|^{2}\approx$
$\displaystyle\left|\left(1-\frac{i\Delta
t}{\hbar}\sum_{j}\langle\psi|\hat{H}_{j}|\psi\rangle-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j}\langle\psi|\hat{H}_{j}^{2}|\psi\rangle-\frac{(\Delta
t)^{2}}{\hbar^{2}}\sum_{j>k}\langle\psi|\hat{H}_{j}\hat{H}_{k}|\psi\rangle\right)\right|^{2}$
(39) $\displaystyle=$ $\displaystyle\left(1-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j}\langle\psi|\hat{H}_{j}^{2}|\psi\rangle-\frac{(\Delta
t)^{2}}{\hbar^{2}}{\tt
Re}\sum_{j>k}\langle\psi|\hat{H}_{j}\hat{H}_{k}|\psi\rangle\right)^{2}$ (40)
$\displaystyle+\left(\frac{\Delta
t}{\hbar}\sum_{j}\langle\psi|\hat{H}_{j}|\psi\rangle+\frac{(\Delta
t)^{2}}{\hbar^{2}}{\tt
Im}\sum_{j>k}\langle\psi|\hat{H}_{j}\hat{H}_{k}|\psi\rangle\right)^{2}.$ (41)

Since

$\displaystyle\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j}\langle\psi|\hat{H}_{j}^{2}|\psi\rangle+\frac{(\Delta
t)^{2}}{\hbar^{2}}{\tt
Re}\sum_{j>k}\langle\psi|\hat{H}_{j}\hat{H}_{k}|\psi\rangle$ (42)
$\displaystyle=$ $\displaystyle\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j}\langle\psi|\hat{H}_{j}^{2}|\psi\rangle+\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sum_{j\neq
k}\langle\psi|\hat{H}_{j}\hat{H}_{k}|\psi\rangle$ (43) $\displaystyle=$
$\displaystyle\frac{(\Delta
t)^{2}}{2\hbar^{2}}\langle\psi|\biggl{(}\sum_{j}\hat{H}_{j}\biggr{)}^{2}|\psi\rangle,$
(44)

we finally have

$\displaystyle|\langle\psi|\phi\rangle|^{2}=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{\hbar^{2}}\langle\psi|\biggl{(}\sum_{j}\hat{H}_{j}\biggr{)}^{2}|\psi\rangle+\frac{(\Delta
t)^{2}}{\hbar^{2}}\biggl{(}\langle\psi|\sum_{j}\hat{H}_{j}|\psi\rangle\biggr{)}^{2}+o((\Delta
t)^{2})$ (45) $\displaystyle=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{\hbar^{2}}\sigma_{\sum\hat{H}_{j}}^{2}(\psi)+o((\Delta t)^{2}).$ (46)

The span of this random walk is thus given by

$\displaystyle\cos s(\phi,\psi)=$
$\displaystyle|\langle\psi|\phi\rangle|+o((\Delta t)^{2})$ (47)
$\displaystyle=$ $\displaystyle 1-\frac{(\Delta
t)^{2}}{2\hbar^{2}}\sigma_{\sum\hat{H}_{j}}^{2}(\psi)+o((\Delta t)^{2}).$ (48)

By using Eqs. (35) and (48), we obtain Eq. (6). Here, we also have the
representation

$\displaystyle\sigma_{\sum\hat{H}_{j}}^{2}(\psi)=$
$\displaystyle\sum_{jk}\left\langle\psi\left|H_{j}H_{k}\right|\psi\right\rangle-\biggl{(}\sum_{j}\left\langle\psi\left|H_{j}\right|\psi\right\rangle\biggr{)}^{2}$
(49) $\displaystyle=$
$\displaystyle\sum_{\chi}\sum_{jk}\left\langle\psi\left|H_{j}|\chi\rangle\langle\chi|H_{k}\right|\psi\right\rangle-\biggl{(}\sum_{j}\left\langle\psi\left|H_{j}\right|\psi\right\rangle\biggr{)}\biggl{(}\sum_{k}\left\langle\psi\left|H_{k}\right|\psi\right\rangle\biggr{)}$
(50) $\displaystyle=$
$\displaystyle\sum_{\chi\neq\psi}\sum_{jk}\left\langle\psi\left|H_{j}|\chi\rangle\langle\chi|H_{k}\right|\psi\right\rangle.$
(51)

##  Appendix C Algorithm for random walk in Hilbert space

A $D$-dimensional quantum state can be parameterized by the $2(D-1)$-real
parameters ${\bm{\theta}}=(\theta_{1},\theta_{2},\cdots,\theta_{2(D-1)})^{\rm
T}$, because of the presence of the overall phase factor and the normalized
condition. The state can be given by

$\displaystyle|\Psi({\bm{\theta}})\rangle=\left(\begin{array}[]{rl}e^{i\theta_{D}}&\cos(\theta_{1})\\\
e^{i\theta_{D+1}}&\sin(\theta_{1})\cos(\theta_{2})\\\
e^{i\theta_{D+2}}&\sin(\theta_{1})\sin(\theta_{2})\cos(\theta_{3})\\\
&\vdots\\\
e^{i\theta_{2(D-1)}}&\sin(\theta_{1})\cdots\sin(\theta_{D-2})\cos(\theta_{D-1})\\\
&\sin(\theta_{1})\cdots\sin(\theta_{D-2})\sin(\theta_{D-1})\end{array}\right).$
(58)

The Fubini-Study metric can be calculated by the relation

$\displaystyle g_{ij}({\bm{\theta}})=$
$\displaystyle\langle\partial_{\theta_{i}}\Psi({\bm{\theta}})|\partial_{\theta_{j}}\Psi({\bm{\theta}})\rangle$
$\displaystyle-\langle\partial_{\theta_{i}}\Psi({\bm{\theta}})|\Psi({\bm{\theta}})\rangle\langle\Psi({\bm{\theta}})|\partial_{\theta_{j}}\Psi({\bm{\theta}})\rangle,$
(59)

for $i,j=1,2,\cdots,2(D-1)$. The distance between states
$|\Psi({\bm{\theta}}+d{\bm{\theta}})\rangle$ and $|\Psi({\bm{\theta}})\rangle$
can be given by

$\displaystyle\Delta
S^{2}=g_{ij}({\bm{\theta}})d\theta_{i}d\theta_{j}=d{\bm{\theta}}^{\rm T}{\bf
g}({\bm{\theta}})d{\bm{\theta}},$ (60)

where $[{\bf g}({\bm{\theta}})]_{ij}=g_{ij}$ and
$d{\bm{\theta}}=(d\theta_{1},d\theta_{2},\cdots,d\theta_{2(D-1)})^{\rm T}$.
Using the form in Eq. (58), we find that $g_{ij}\in{\mathbb{R}}$ for $i,j\leq
D-1$ and $i,j\geq D$. We also have $(g_{ij}+g_{ji})d\theta_{i}d\theta_{j}=0$.
We can thus take ${\bf g}$ as a real symmetric matrix by replacing
$g_{ij},g_{ji}$ with $(g_{ij}+g_{ji})/2$.

We then consider the uniform random walk in the Hilbert space, where each step
length is fixed with the Fubini-Study distance $\Delta S$. We first generate
an initial state $|\Psi_{t=0}\rangle$, which is constracted with the random
phase and the random amplitude that satisfies the normalization condition. The
parameters $\theta_{i=1,\cdots,2(D-1)}$ are obtained by inversely solving Eq.
(58). Since the parameters ${\bm{\theta}}$ are given in the non-orthogonal
coordinates with the curved space, we consider the local coordinate
transformation to the local orthogonal coordinate at the point ${\bm{\theta}}$
in order to consider the uniformly distributed random walk in the Hilbert
space. In this case, the distance $\Delta S$ can be given by

$\displaystyle\Delta S^{2}=$ $\displaystyle\Delta{\bm{\Theta}}^{\rm T}{\bf
h}\Delta{\bm{\Theta}}=\sum_{i=1}^{2(D-1)}h_{i}\Delta\Theta_{i}^{2},$ (61)

where we used the fact that the real symmetric matrix ${\bf g}$ is
diagonalized by orthogonal matrix ${\bf V}$ as ${\bf V}^{\rm T}{\bf g}{\bf
V}={\bf h}={\rm diag}[h_{1},h_{2},\cdots,h_{2(D-1)}]$ with
$h_{i}\in{\mathbb{R}}$, ${\bf g}{\bf v}_{i}=h_{i}{\bf v}_{i}$, ${\bf V}=({\bf
v}_{1},{\bf v}_{2},\cdots,{\bf v}_{2(D-1)})$, and $\Delta{\bm{\Theta}}={\bf
V}^{\rm T}\Delta{\bm{\theta}}$. Since ${\bf v}_{i}\in{\mathbb{R}}^{2(D-1)}$,
we also have $\Delta{\bm{\Theta}}\in{\mathbb{R}}^{2(D-1)}$. We generate the
vector ${\bf u}=(u_{1},u_{2},\cdots,u_{2(D-1)})^{\rm T}$ uniformly distributed
over the surface of the unit $2(D-1)$-ball, which can easily constructed from
the $2(D-1)$-real random values $(x_{1},x_{2},\cdots,x_{2(D-1)})$ generated
from the normal distribution ${\mathcal{N}}(0,1)$ with the relation
$u_{i}=x_{i}/\sqrt{x_{1}^{2}+x_{2}^{2}+\cdots+x_{2(D-1)}^{2}}$. Using the
uniformly distributed unit vector ${\bf u}$, we can take

$\displaystyle\Delta\Theta_{i}=\frac{\Delta S}{\sqrt{h_{i}}}u_{i},$ (62)

which satisfies the relation (61). The inverse transformation provides
$\Delta{\bm{\theta}}={\bf V}\Delta{\bm{\Theta}}$.

By employing a state $|\Psi_{t}({\bm{\theta}})\rangle$, we can generate a new
state $|\Psi_{t+1}({\bm{\theta}}+\Delta{\bm{\theta}})\rangle$ from one of the
uniformly distributed random vectors ${\bf u}$, the distance of which between
$|\Psi_{t}\rangle$ and $|\Psi_{t+1}\rangle$ is $\Delta S$. The Fubini-Study
distance between the state $|\Psi_{t}\rangle$ and the initial state
$|\Psi_{t=0}\rangle=|0\rangle$ is given by
$L_{t}=\cos^{-1}(|\langle\Psi_{t}|0\rangle|)$. At the random walk step $t=M$,
if the state $|\Psi_{t=M}\rangle$ reaches the distance
$\pi/2-L_{t=M}\leq\epsilon$ with a small parameter, we consider that the state
starting from $|\Psi_{t=0}\rangle=|0\rangle$ arrives a farthest state through
the random walk.

##  Appendix D Concentration Inequality

Let $|\Psi\rangle=(\Psi_{1},\cdots,\Psi_{D})^{\rm T}$ and
$|\phi\rangle=(\phi_{1},\cdots,\phi_{D})^{\rm T}$ be $D$-dimensional random
normalized complex vectors. Here, $\Psi_{i}$ and $\phi_{i}$ for $i=1,\cdots,D$
are generated from mutually independent random variables
$u_{j},v_{j}\in{\mathbb{R}}$ for $j=1,\cdots 2D$, where
$\Psi_{j}=(u_{2j-1}+iu_{2j})$ and $\phi_{j}=(v_{2j-1}+iv_{2j})$ with
$\sum\limits_{j=1}^{2D}u_{j}^{2}=\sum\limits_{j=1}^{2D}v_{j}^{2}=1$. The
fidelity between two states is given by $|\langle\phi|\Psi\rangle|^{2}=({\bf
v}\cdot{\bf u})^{2}+({\bf v}^{\prime}\cdot{\bf u})^{2}$, where ${\bf
u}=(u_{1},u_{2},\cdots,u_{2D})^{\rm T}$, ${\bf
v}=(v_{1},v_{2},\cdots,v_{2D})^{\rm T}$ and ${\bf
v}^{\prime}=(-v_{2},v_{1},-v_{4},v_{3},\cdots,-v_{2D},v_{2D-1})^{\rm T}$. If
${\bf u}$ and ${\bf v}$ are $2D$-dimensional isotropically random unit
vectors, the angular distance is known to show the relation $\cos^{-1}(|{\bf
v}\cdot{\bf u}|)=\pi/2+{\mathcal{O}}(1/\sqrt{2D})$ [21]. As a result, in
integrated qubit systems, two random states generated isotropically are almost
orthogonal. We naturally expect that the fidelity between two random states
generated isotropically is exponentially small with respect to the number of
qubits $N$:

$\displaystyle|\langle\phi|\Psi\rangle|^{2}={\mathcal{O}}\left(\frac{1}{2^{N}}\right).$
(63)

For an interesting example, we consider the concentration inequality for the
fidelity between a basis $|i\rangle$ and an isotropically random states
$|\phi\rangle=\sum\limits_{i=1}^{D}\phi_{i}|i\rangle$ in a $D$-dimensional
Hilbert space with $\phi_{j}=(v_{2j-1}+iv_{2j})$ for ${\bf
v}=(v_{1},\cdots,v_{2D})^{\rm T}\in{\mathbb{R}}^{2D}$. The vector ${\bf v}$ is
a $2D$-dimensional isotropically random vector, the components of which are
mutually independent and can be generated from the $2D$-dimensional normal
distribution ${\bm{\xi}}=(\xi_{1},\cdots,\xi_{2D})^{\rm
T}\sim{\mathcal{N}}_{2D}({\bf 0},{\bf 1})$ with a relation
$v_{i}=\xi_{i}/|{\bm{\xi}}|$ and
$|{\bm{\xi}}|^{2}=\sum\limits_{i}^{2D}\xi_{i}^{2}$. In this case, the fidelity
is given by

$\displaystyle|\langle
i|\phi\rangle|^{2}=v_{2i-1}^{2}+v_{2i}^{2}=\frac{\xi_{2i-1}^{2}+\xi_{2i}^{2}}{|{\bm{\xi}}|^{2}}.$
(64)

We here consider the following probability

$\displaystyle P\left(\left||\langle
i|\phi\rangle|^{2}-\frac{1}{D}\right|\geq\epsilon\right).$ (65)

This probability can be given by

$\displaystyle P\left(\left||\langle
i|\phi\rangle|^{2}-\frac{1}{D}\right|\geq\epsilon\right)$ (66)
$\displaystyle\leq$ $\displaystyle P\left(|\langle i|\phi\rangle|^{2}\geq
A_{+}\right)+P\left(|\langle i|\phi\rangle|^{2}\leq A_{-}\right),$ (67)

where $A_{\pm}=1/D\pm\epsilon$.

By using the Fréchet inequality, the first term can be reduced to

$\displaystyle P_{+}\equiv$ $\displaystyle P\left(|\langle
i|\phi\rangle|^{2}\geq A_{+}\right)$ (68) $\displaystyle=$ $\displaystyle
P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\geq|{\bm{\xi}}|^{2}A_{+})$ (69)
$\displaystyle\leq$ $\displaystyle{\rm
min}\\{P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\geq
2DA_{+}\eta_{+}),P(2D\eta_{+}\geq|{\bm{\xi}}|^{2})\\},$ (70)

for $\eta_{+}>0$.

The first part can be given by

$\displaystyle P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\geq 2DA_{+}\eta_{+})$ (71)
$\displaystyle=$ $\displaystyle
P\left(\frac{\xi_{2i-1}^{2}+\xi_{2i}^{2}}{2}-1\geq DA_{+}\eta_{+}-1\right)$
(72) $\displaystyle\leq$ $\displaystyle
P\left(\left|\frac{\xi_{2i-1}^{2}+\xi_{2i}^{2}}{2}-1\right|\geq
DA_{+}\eta_{+}-1\right).$ (73)

It is known that if $Z_{k}$ is generated from the normal distribution
${\mathcal{N}}(0,1)$, the concentration of $\chi^{2}$-variables holds:

$\displaystyle P\left(\left|\frac{1}{n}\sum_{k=1}^{n}Z_{k}^{2}-1\right|\geq
t\right)\leq 2\exp\left(-\frac{nt^{2}}{8}\right),$ (74)

for $t\in(0,1)$ [22]. Using this concentration inequality, we have a tail
bound:

$\displaystyle P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\geq 2DA_{+}\eta_{+})\leq$
$\displaystyle 2\exp\left[-\frac{(DA_{+}\eta_{+}-1)^{2}}{4}\right],$ (75)

for $1/(1+D\epsilon)<\eta_{+}<2/(1+D\epsilon)$.

By using the tail bound (74), the second part can be given by

$\displaystyle P(2D\eta_{+}\geq|{\bm{\xi}}|^{2})=$ $\displaystyle
P(2D-|{\bm{\xi}}|^{2}\geq 2D(1-\eta_{+}))$ (76) $\displaystyle\leq$
$\displaystyle P\left(\left|\frac{|{\bm{\xi}}|^{2}}{2D}-1\right|\geq
1-\eta_{+}\right)$ (77) $\displaystyle\leq$ $\displaystyle
2\exp\left[-\frac{D(1-\eta_{+})^{2}}{4}\right],$ (78)

where $0<\eta_{+}<1$.

As a result, we have

$\displaystyle P_{+}\leq$ $\displaystyle{\rm
min}\biggl{\\{}\inf\limits_{\eta_{+}\in
R_{+}}2\exp\left(-\frac{X_{+}}{4}\right),\inf\limits_{\eta_{+}\in
R_{+}}2\exp\left(-\frac{Y_{+}}{4}\right)\biggr{\\}},$ (79)

where $X_{+}\equiv[(1+D\epsilon)\eta_{+}-1]^{2}$, $Y_{+}\equiv
D(1-\eta_{+})^{2}$, and

$\displaystyle R_{+}:\frac{1}{1+D\epsilon}<\eta_{+}<{\rm
min}\left\\{1,\frac{2}{1+D\epsilon}\right\\}.$ (80)

A simple computation shows

$\displaystyle\inf\limits_{\eta_{+}\in
R_{+}}2\exp\left(-\frac{X_{+}}{4}\right)=$ $\displaystyle
2\exp\left[-\frac{1}{4}{\rm min}\\{1,(\epsilon D)^{2}\\}\right],$
$\displaystyle\inf\limits_{\eta_{+}\in
R_{+}}2\exp\left(-\frac{Y_{+}}{4}\right)=$ $\displaystyle
2\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right],$

and it is simply shown that the probability $P_{+}$ is bounded as

$\displaystyle P_{+}\leq
2\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right].$
(81)

We now consider the second term $P_{-}\equiv P\left(|\langle
i|\phi\rangle|^{2}\leq A_{-}\right)$. For $1/D<\epsilon$, we have
$P\left(|\langle i|\phi\rangle|^{2}\leq A_{-}<0\right)=0$. For $\epsilon<1/D$.
we consider

$\displaystyle P_{-}=$ $\displaystyle
P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\leq|{\bm{\xi}}|^{2}A_{-})$ (82)
$\displaystyle\leq$ $\displaystyle{\rm
min}\\{P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\leq
2D\eta_{-}A_{-}),P(2D\eta_{-}\leq|{\bm{\xi}}|^{2})\\},$ (83)

for $\eta_{-}>0$. The first part can be given by

$\displaystyle P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\leq 2D\eta_{-}A_{-})$ (84)
$\displaystyle=$ $\displaystyle
P\left(\frac{\xi_{2i-1}^{2}+\xi_{2i}^{2}}{2}-1\leq D\eta_{-}A_{-}-1\right)$
(85) $\displaystyle\leq$ $\displaystyle
P\left(\left|\frac{\xi_{2i-1}^{2}+\xi_{2i}^{2}}{2}-1\right|\geq
1-DA_{-}\eta_{-}\right).$ (86)

Using concentration inequality (74), we have:

$\displaystyle P(\xi_{2i-1}^{2}+\xi_{2i}^{2}\leq 2DA_{-}\eta_{-})\leq$
$\displaystyle 2\exp\left[-\frac{(1-DA_{-}\eta_{-})^{2}}{4}\right],$ (87)

for $0<\eta_{-}<1/(1-D\epsilon)$. By using the two-sided tail bound (74), the
second part can be given by

$\displaystyle P(2D\eta_{-}\leq|{\bm{\xi}}|^{2})\leq$ $\displaystyle
P\left(\left|\frac{1}{2D}|{\bm{\xi}}|^{2}-1\right|\geq\eta_{-}-1\right)$ (88)
$\displaystyle\leq$ $\displaystyle
2\exp\left[-\frac{D(\eta_{-}-1)^{2}}{4}\right],$ (89)

where $1<\eta_{-}<2$.

As a result, we have

$\displaystyle P_{-}\leq$ $\displaystyle{\rm
min}\biggl{\\{}\inf\limits_{\eta_{-}\in
R_{-}}2\exp\left(-\frac{X_{-}}{4}\right),\inf\limits_{\eta_{-}\in
R_{-}}2\exp\left(-\frac{Y_{-}}{4}\right)\biggr{\\}},$ (90)

where $X_{-}\equiv[1-\eta_{-}(1-D\epsilon)]^{2}$, $Y_{-}\equiv
D(\eta_{-}-1)^{2}$, and

$\displaystyle R_{-}:1<\eta_{-}<{\rm
min}\left\\{2,\frac{1}{1-D\epsilon}\right\\}.$ (91)

A simple computation shows

$\displaystyle\inf\limits_{\eta_{-}\in
R_{-}}2\exp\left(-\frac{X_{-}}{4}\right)=$ $\displaystyle
2\exp\left[-\frac{(\epsilon D)^{2}}{4}\right],$
$\displaystyle\inf\limits_{\eta_{-}\in
R_{-}}2\exp\left(-\frac{Y_{-}}{4}\right)=$ $\displaystyle
2\exp\left[-\frac{D}{4}{\rm
min}\left\\{1,\left(\frac{D\epsilon}{1-D\epsilon}\right)^{2}\right\\}\right]$
$\displaystyle\leq$ $\displaystyle
2\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right],$

and it is simply shown that the probability $P_{-}$ is also bounded as

$\displaystyle P_{-}\leq
2\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right].$
(92)

Using (67), (81) and (92), we finally obtain the concentration inequality

$\displaystyle P\left(\left||\langle
i|\phi\rangle|^{2}-\frac{1}{D}\right|\geq\epsilon\right)\leq
4\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right].$
(93)

We here assume that an arbitrary complex vector $|\Psi\rangle$ is generated
from a basis $|i\rangle$ by using a unitary matrix $\hat{U}$, where
$|\Psi\rangle=\hat{U}|i\rangle$. Using the same unitary matrix, the state
$|\phi\rangle$ is transformed to
$|\Psi^{\prime}\rangle=\hat{U}|\phi\rangle=(\Psi_{1}^{\prime},\Psi_{2}^{\prime},\cdots,\Psi_{D}^{\prime})^{\rm
T}$, with $\Psi_{i}^{\prime}=(\zeta_{2i-1}+i\zeta_{2i})/|{\bm{\zeta}}|$ where
$|{\bm{\zeta}}|=\sum_{i=1}^{2D}\zeta_{i}^{2}$. If we write
$|\phi\rangle=|{\bm{\xi}}\rangle/\sqrt{\langle{\bm{\xi}}|{\bm{\xi}}\rangle}$
and
$|\Psi^{\prime}\rangle=|{\bm{\zeta}}\rangle/\sqrt{\langle{\bm{\zeta}}|{\bm{\zeta}}\rangle}$,
we have $|{\bm{\zeta}}\rangle=\hat{U}|{\bm{\xi}}\rangle$ and
$\sum\limits_{i=1}^{2D}\xi_{i}^{2}=\langle{\bm{\xi}}|{\bm{\xi}}\rangle=\langle{\bm{\zeta}}|{\bm{\zeta}}\rangle=\sum\limits_{i=1}^{2D}\zeta_{i}^{2}$,
due to the unitary transformation. In this case, if
${\bm{\xi}}=(\xi_{1},\cdots,\xi_{2D})^{\rm T}$ is generated from the
$2D$-dimensional normal distribution ${\bm{\xi}}\sim{\mathcal{N}}_{2D}({\bf
0},{\bf 1})$, the $2D$-dimensional random vector
${\bm{\zeta}}=(\zeta_{1},\cdots,\zeta_{2D})^{\rm T}$ is also uniformly
generated from the $2D$-dimensional normal distribution
${\bm{\zeta}}\sim{\mathcal{N}}_{2D}({\bf 0},{\bf 1})$.

To conclude, an arbitrary state $|\Psi\rangle$ and an isotropically random
state $|\Psi^{\prime}\rangle$ is almost orthogonal for a large dimension
Hilbert space, the concentration inequality of which is given by

$\displaystyle
P\left(\left||\langle\Psi|\Psi^{\prime}\rangle|^{2}-\frac{1}{D}\right|\geq\epsilon\right)\leq
4\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right].$
(94)

or

$\displaystyle
1-4\exp\left[-\frac{D}{4}\left(\frac{D\epsilon}{1+D\epsilon}\right)^{2}\right]\leq
P\left(\left||\langle\Psi|\Psi^{\prime}\rangle|^{2}-\frac{1}{D}\right|<\epsilon\right)\leq
1.$ (95)

##  Appendix E Percolation model: Numerical approach

We first randomly prepare $M$-complex vectors
$|\Psi_{n}\rangle=(u_{1}^{(n)},u_{2}^{(n)},\cdots,u_{D}^{(n)})^{\rm T}$ for
$n=1,2,\cdots M$, where $D$ is the dimension of the Hilbert space related to
the number of qubits $N$ as $D=2^{N}$. For uniformly preparing normalized
states at random in the Hilbert space, we generate $2D$-real values
$(x_{1},x_{2},\cdots,x_{2D})$ from the normal distribution
${\mathcal{N}}(0,1)$, and provide $u_{j}^{(n)}=(x_{2j-1}+ix_{2j})/r$ for
$j=1,2,\cdots,D$, where $r=\sqrt{x_{1}^{2}+x_{2}^{2}+\cdots+x_{2D}^{2}}$. We
then evaluated the Fubini–Study distance
$d_{mn}=\cos^{-1}(|\langle\Psi_{m}|\Psi_{n}\rangle|)$ between all pairs of
random states for $m,n=1,2,\cdots,M$, where the Fubini-Study distance ranges
from $0$ to $\pi/2$.

By using the following procedures, we construct clusters of states that are
linked with the condition where the Fubini–Study distance $d_{mn}$ is less
than a given $\Delta S$. After initially preparing boolean array $b_{mn}={\tt
F}$ for all pairs of states, we replace $b_{mn}={\tt F}$ with ${\tt T}$ if
$d_{mn}\leq dS$ for $m,n=1,2,\cdots,M$. We then make $b_{ml},b_{nl}=b_{ml}\lor
b_{nl}$ for $l=1,\cdots,M$, if $b_{mn}={\tt T}$, where $b_{ml}$ and $b_{nl}$
share the information of states which are linked from the states $m$ and $n$.
The boolean array $b_{mn}$ has then all information of random states that
belongs to a cluster; if we find $b_{mn}={\tt T}$ for $n(\neq m)$, a state
$|\Psi_{n}\rangle$ belongs to a cluster with a state $|\Psi_{m}\rangle$. Let
$N_{\rm c}$ be the number of clusters generated by the rule shown above. In
order to minimize the duplication of cluster information, if $b_{ml}=b_{nl}$
for $m>n$, we replace $b_{ml}$ with ${\tt F}$. We then generate the boolean
array for the cluster $c_{\alpha}=\\{b_{1n},b_{2n},\cdots,b_{Mn}\\}$, if at
least one $b_{mn}$ is ${\tt T}$ for $m=1,\cdots,M$. The index $\alpha$
indicates the label of the cluster, where $\alpha$ corresonds to one-to-one
with $n$ such that $b_{mn}={\tt T}$ for ${}^{\exists}m$. If
$b_{ln}=b_{mn}={\rm T}$ for $b_{ln,mn}\in c_{\alpha}$, the states
$|\Psi_{l}\rangle$ and $|\Psi_{m}\rangle$ belong to the $\alpha$-th cluster,
the distance of which is $d_{lm}$. By comparing all the distances between
states in the $\alpha$-th cluster, we can find the maximum value of $d_{lm}$.
We call this value the maximum distance $L_{\alpha}$ in the $\alpha$-the
cluster. If we find $\pi/2-L_{\alpha}\leq\epsilon$ with a small parameter
$\epsilon$, we regard that the $\alpha$-th cluster spans maximally in the
Hilbert space.