Finding spin glass ground states using quantum walks

A continuous-time quantum walk is defined on an undirected graph G(V, E), with $V=\{j\}{}_{j=0}^{N-1}$ the set of N vertex labels and E the set of label-pairs (j, k) associated with edges. The vertices correspond to the positions of the walker, and the edges indicate the allowed transitions between vertices. This is conveniently encoded in the adjacency matrix A of the graph, which has entries A_jk = 1 for $(j,k)\in E$ and A_jk = 0 otherwise. The Laplacian of G is $L=A-D$, where D is a diagonal matrix formed from the degree of each vertex, ${D}_{{jj}}={\rm{\deg }}(j)$, where ${\rm{\deg }}(j)$ is the number of edges connected to vertex j. Both the adjacency matrix A and Laplacian L are symmetric matrices which can thus be used to define a quantum Hamiltonian for the dynamics of the continuous-time quantum walk on the graph. In this work, we only need regular graphs, for which ${\rm{\deg }}(j)$ is constant with respect to j. For regular graphs, the only difference between using the adjacency matrix A or Laplacian L is an irrelevant global phase (Childs and Goldstone 2004). We use the Laplacian form of the Hamiltonian for consistency with prior work. We thus define the quantum walk Hamiltonian ${\hat{H}}_{G}$ for a quantum walk on graph G by

$$\begin{eqnarray}&&\left\langle j\right|{\hat{H}}_{G}\left|k\right\rangle =-\gamma {L}_{{jk}},\end{eqnarray} \tag{ 5 }$$

where γ is the hopping rate between connected vertices per unit time. The states $\left|j\right\rangle ,\left|k\right\rangle$ for $j,k\in V$ are associated with the vertices of G and form a basis for a Hilbert space of dimension N. In the Ising model context, the dimension of the Hilbert space is $N={2}^{n}$ where n is the number of qubits, and $\{| j\rangle \}{}_{j=0}^{N-1}$ is the computational basis. For a quantum walk starting in state $\left|\psi (0)\right\rangle$, the state of the walker evolves according to the Schrödinger equation, with formal solution

$$\begin{eqnarray}&&\left|\psi (t)\right\rangle =\exp \{-{\rm{i}}{\hat{H}}_{G}t\}\left|\psi (0)\right\rangle ,\end{eqnarray} \tag{ 6 }$$

using units in which ${\hslash }=1$.

The task is to solve an optimization problem whose N = 2ⁿ candidate solutions j are represented in the computational basis $\{\left|j\right\rangle \}{}_{j=0}^{N-1}$, where j is a bit string corresponding to the state of n qubits. The problem is encoded in an Ising Hamiltonian ${\hat{H}}_{P}$, of the form described by ${\hat{H}}_{I}$ in (1) and whose eigenbasis is the computational basis. We write the basis state with eigenvalue ${E}_{a}^{\left(P\right)}$ as $\left|{E}_{a}^{\left(P\right)}\right\rangle$, with $a\in \{0...N-1\}$, and adopt the convention that ${E}_{a}^{\left(P\right)}\leqslant {E}_{a+1}^{\left(P\right)}$. In other words, ${\left\{\left|{E}_{a}^{\left(P\right)}\right\rangle \right\}}_{a=0}^{N-1}$ is a reordering of $\{\left|j\right\rangle \}{}_{j=0}^{N-1}$ based on the corresponding eigenenergies of ${\hat{H}}_{P}$. The encoding is chosen such that the solution corresponds to the ground state $\left|{E}_{0}^{\left(P\right)}\right\rangle$ of the problem Hamiltonian ${\hat{H}}_{P}$.

To use a quantum walk to solve the problem, we must first choose a suitable state in which to initialize the system. With no prior knowledge of the solution, the equal superposition of all basis states

$$\begin{eqnarray}&&\left|\psi (0)\right\rangle ={N}^{-1/2}\displaystyle \sum _{j=0}^{N-1}\left|j\right\rangle ,\end{eqnarray} \tag{ 7 }$$

is a sensible choice that avoids bias. More generally, the initial state can be prepared as weighted or biased superposition, to incorporate prior knowledge about the solution (Perdomo-Ortiz et al 2011, Duan et al 2013, Chancellor 2017, Graß and Lewenstein 2017, Baldwin and Laumann 2018, Kechedzhi et al 2018, Graß 2019). Next, we choose a suitable walk graph G. The main requirement is that the ground state of the quantum walk Hamiltonian ${\hat{H}}_{G}$ coincides with the initial state, either biased or unbiased (see section 6.2). A simple way to achieve a biased starting state would be to 'tilt' the driver fields so they are no longer completely transverse. We only treat the unbiased case in this work, so our initial state will be $\left|\psi (0)\right\rangle$ throughout. The full Hamiltonian $\hat{H}(\gamma )$ is defined by adding the quantum walk Hamiltonian ${\hat{H}}_{G}$ to the problem Hamiltonian ${\hat{H}}_{P}$

$$\begin{eqnarray}&&\hat{H}(\gamma )\equiv {\hat{H}}_{G}+{\hat{H}}_{P},\end{eqnarray} \tag{ 8 }$$

where the key parameter is the hopping rate γ in ${\hat{H}}_{G}$, see (5). The computation is performed by evolving the initial state (7) under the full Hamiltonian $\hat{H}(\gamma )$ for a time t_f, then measuring the qubit register in the computational basis. The intuition, based on the faster spreading of quantum walks over classical found in prior work (Farhi and Gutmann 1998), is that the quantum walk dynamics provide rapid exploration of the basis states, while the energy structure of the problem Hamiltonian ${\hat{H}}_{P}$ causes localisation around low-energy states.

The success probability $P({t}_{f})={\left|\left\langle {E}_{0}^{\left(P\right)}| \psi ({t}_{f})\right\rangle \right|}^{2}$ of finding the solution state when measuring will not in general be unity. It will typically be necessary to repeat the protocol multiple times to obtain a high probability of success over all the repeats. In general, it will be best to use different measurement times t_f for each repeat. Different measurement times will produce different success probabilities P(t_f), and varying the measurement time avoids repeatedly measuring at a time for which the probability P(t_f) happens to be atypically small. More precisely, we choose the measurement time t_f uniformly at random in an interval $\left[t,t+{\rm{\Delta }}t\right]$, and define an average single run success probability

$$\begin{eqnarray}&&\bar{P}(t,{\rm{\Delta }}t)\equiv \displaystyle \frac{1}{{\rm{\Delta }}t}{\int }_{t}^{t+{\rm{\Delta }}t}{\rm{d}}{t}_{f}P({t}_{f}).\end{eqnarray} \tag{ 9 }$$

Operationally, choosing the measurement time t_f randomly in the interval $[t,t+{\rm{\Delta }}t]$ samples success probabilities from the distribution with $\bar{P}(t,{\rm{\Delta }}t)$ as its mean. Sampling measurement times in this way means that the protocol typically needs to be repeated ${M}_{\mathrm{rep}}\sim 1/\bar{P}(t,{\rm{\Delta }}t)$ times to achieve an overall O(1) success probability. Note that it is not generally possible to check whether the state measured is indeed the ground state of ${\hat{H}}_{P}$. However, it is easy to calculate the energy of the state measured in each repeat. If only the lowest energy state is accepted, it is only necessary for the ground state of ${\hat{H}}_{P}$ to be measured once out of all the repeats. The more repeats, the more confidence is gained that the lowest energy state found is the ground state. And studying the distribution of the sampled energies can provide more information about the problem.

The procedure described in this subsection does not in general provide an optimal quantum algorithm, because the repeats do not use information gained from the outcomes of previous runs. We will discuss this further in section 7; for most of this paper we are concerned with understanding the average single run success probability, as an essential prerequisite to building optimal algorithms.

In the limit of small interval width Δt, the average success probability defined in (9) reduces to the single time probability $P({t}_{f})={\mathrm{lim}}_{{\rm{\Delta }}t\to 0}\bar{P}({t}_{f},{\rm{\Delta }}t)$. The long time limit of this average,

$$\begin{eqnarray}&&{P}_{\infty }\equiv \bar{P}(0,\infty )\equiv \mathop{\mathrm{lim}}\limits_{{\rm{\Delta }}t\to \infty }\bar{P}(0,{\rm{\Delta }}t),\end{eqnarray} \tag{ 10 }$$

is particularly useful, because it can be calculated via a numerical diagonalization of the Hamiltonian (see section 4) and it predicts the short time average well (see section 5.3). In this paper, we will often use the long time average ${P}_{\infty }$ as an indication of the success probability achievable in a single run, and thus the number of repeats required to achieve O(1) success probability overall. We will separately address the timescale required to reach this probability in each run.

There are many graph-based Hamiltonians with the initial state $\left|\psi (0)\right\rangle$ defined in (7) as the ground state. A common choice is the complete graph K, in which every vertex is connected to every other. This graph has the quantum walk Hamiltonian ${\hat{H}}_{K}$ that couples every computational basis state $\left|j\right\rangle$ state to every other,

$$\begin{eqnarray}\begin{array}{rcl}{\hat{H}}_{K} & = & \gamma \left[N{\mathbb{1}}-\displaystyle \sum _{j,k=0}^{N-1}\left|k\right\rangle \left\langle j\right|\right]\\ & = & \gamma N\left[{\mathbb{1}}-\left|\psi (0)\right\rangle \left\langle \psi (0)\right|\right].\end{array}\end{eqnarray} \tag{ 11 }$$

The complete graph is useful because it makes some algorithms analytically tractable (see, e.g. Childs and Goldstone 2004). However, for implementation on qubit-based hardware, the complete graph is not in general practical, requiring higher order interaction terms than the transverse Ising term (3). In this qubit setting, an implementation of the complete graph requires a sum over every one-body term (e.g. ${\hat{X}}_{j}$), every two-body term (e.g. ${\hat{X}}_{j}{\hat{X}}_{k}$), every three-body term (e.g. ${\hat{X}}_{j}{\hat{X}}_{k}{\hat{X}}_{l}$) ... up to the n-body term ${\prod }_{j=0}^{n-1}{\hat{X}}_{j}$, a total of N terms. One- and two-body terms are relatively easy to implement, since they correspond to Hamiltonians found naturally. Terms in three or more Pauli-X operators are much more difficult and generally require extra qubits to engineer in real physical systems.

A more natural choice of graph for qubits is the hypercube. The n-bit labels are associated with the vertices of the graph such that the edges correspond to flipping one bit, as illustrated in figure 1. The hypercube quantum walk Hamiltonian ${\hat{H}}_{h}$ on n qubits is composed of single-body terms

$$\begin{eqnarray}&&{\hat{H}}_{h}=\gamma \left[n{\mathbb{1}}-\displaystyle \sum _{j=0}^{n-1}{\hat{X}}_{j}\right].\end{eqnarray} \tag{ 12 }$$

With ${\hat{H}}_{h}$ as the graph Hamiltonian, the full quantum walk computational Hamiltonian $\hat{H}(\gamma )$ defined in (8) is a transverse Ising Hamiltonian in the form of ${\hat{H}}_{{TI}}$ in (4), with the control functions A(t) and B(t) kept constant throughout the computation.

In this work, we predominantly use the hypercube graph, with some comparisons made with the same problems on the complete graph.

The simplest example of an algorithm in this continuous-time quantum walk setting is the search problem. The problem is to find the marked state, a single bit-string $m\in \{0,1\}{}^{n}$ out of N = 2ⁿ possible bit strings. Finding a marked state was shown to have a quantum algorithm with a speed up over classical algorithms by Grover (1996). To map this problem to the continuous-time Hamiltonian setting, the marked basis state $\left|m\right\rangle$ is given one less unit of energy than all the rest of the basis states, by defining the problem Hamiltonian ${\hat{H}}_{{\rm{S}}}$ as

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{S}}}=-| m\rangle \langle m| .\end{eqnarray} \tag{ 13 }$$

By construction, the problem Hamiltonian ${\hat{H}}_{{\rm{S}}}$ has the marked state $| m\rangle$ as its ground state.

The continuous-time quantum walk search problem has been analytically solved (Childs and Goldstone 2004) for several different walk graphs. For the complete graph and the hypercube graph, a quantum speed up is obtained for carefully chosen optimal values of the hopping rate γ. For the complete graph Hamiltonian ${\hat{H}}_{K}$, the optimal value is ${\gamma }_{\mathrm{opt}}^{(K)}=1/N$, while for the hypercube Hamiltonian, ${\hat{H}}_{h}$, the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ is given by

$$\begin{eqnarray}&&2{\gamma }_{\mathrm{opt}}^{(h)}=\displaystyle \frac{1}{N}\displaystyle \sum _{r=1}^{n}\displaystyle \left(\genfrac{}{}{0em}{}{n}{r}\right)\displaystyle \frac{1}{r},\end{eqnarray} \tag{ 14 }$$

where $\displaystyle \left(\genfrac{}{}{0em}{}{n}{r}\right)=\tfrac{n!}{r!(n-r!)}$ is the binomial coefficient. For a quantum speed up, the hopping rate must be set to ${\gamma }_{\mathrm{opt}}^{(h)}$ as defined by (14) with high precision. It has been shown (Morley 2019) that the fractional tolerance to misspecification of the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ falls as $O({N}^{-1/2})$.

The measurement time must also be chosen appropriately. In the limit of large problem size N, the marked state can be found with unit success probability, ${\mathrm{lim}}_{N\to \infty }\left[P({t}_{f}^{(\mathrm{opt})})\right]=1$, by measuring in the computational basis at an optimal measurement time ${t}_{f}^{(\mathrm{opt})}$. For both the hypercube and complete graphs, the optimal time ${t}_{f}^{(\mathrm{opt})}$ scales with the square-root of the problem size N as ${t}_{f}^{(\mathrm{opt})}\simeq \tfrac{\pi }{2}{N}^{1/2}$. This corresponds to a quadratic speed up compared to the best classical algorithm. Due to the absence of structure in the search problem specifically, such a quadratic speed up has been proven to the best possible quantum speed up (Bennett et al 1997).

The variation of P(t_f) with t_f is shown in figure 2(a) for search on hypercube graphs of size N = 2³⁰ (i.e. n = 30 qubits) and N = 2¹¹ (i.e. n = 11 qubits), using the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$. The sinusoidal oscillations of the probability P(t_f) occur because the quantum walk is performing Rabi oscillations between the initial state and the marked state. The two lowest energy levels of the full Hamiltonian $\hat{H}(\gamma )$ with varying γ undergo an avoided level crossing at ${\gamma }_{\mathrm{opt}}^{(h)}$ and the associated eigenstates $\left|{E}_{0}({\gamma }_{\mathrm{opt}}^{(h)}\right\rangle )$ and $\left|{E}_{1}({\gamma }_{\mathrm{opt}}^{(h)}\right\rangle )$ are approximately the orthogonal equal superpositions of the starting state and marked state, $\left|{E}_{\mathrm{0,1}}({\gamma }_{\mathrm{opt}}^{(h)}\right\rangle )\,\simeq (\left|\psi (0)\right\rangle \pm \left|m\right\rangle )/{2}^{\tfrac{1}{2}}$. The gap ${E}_{1}({\gamma }_{\mathrm{opt}}^{(h)})-{E}_{0}({\gamma }_{\mathrm{opt}}^{(h)})$ scales with the problem size N as $O({N}^{-1/2})$ (Childs and Goldstone 2004).

Zoom In Zoom Out Reset image size

These simple, two-level dynamics describe the quantum walk solution to the search problem well for large problem size N: the oscillations in the N = 2³⁰ case have no visible irregularities. For smaller sizes, finite-size effects due to population of higher energy levels are apparent: the oscillations in the $N={2}^{11}$ case have lower probability peaks and show some irregular behaviour, such as the small dip on the first peak. These finite-size effects are further illustrated in figure 2(b), which shows the instantaneous success probability P(t_f) at the asymptotically optimal and numerically determined best times, as well as the infinite-time average success probability ${P}_{\infty }$ defined in (10). All three probabilities show a pronounced dip around n = 8 qubits, with smooth behaviour only settling in for n > 12 qubits. Figure 2(b) also shows that the infinite-time probability ${P}_{\infty }$ asymptotes to a half. Hence, a quantum walk search with a random measurement time should on average only need to be repeated twice to locate the marked state; knowing the exact time to measure for the optimal success probability is not necessary for the success of the algorithm. Fixed point quantum search algorithms (Yoder et al 2014, Dalzell et al 2017) are another approach that avoids the need to know how long to run the algorithm for.

The search problem in the continuous-time quantum computing setting has two important drawbacks. Firstly, implementing the problem Hamiltonian ${\hat{H}}_{S}$ directly on n qubits requires O(2ⁿ) terms of products of up to n Pauli-Z operators, similar to the problem with implementing the complete-graph Hamiltonian ${\hat{H}}_{K}$, defined in (11), on qubits. Implementing higher order Pauli-Z terms can be done using extra qubits as 'gadgets', e.g. Jordan and Farhi (2008). An alternative type of gadget, specifically for permutation-symmetric problems like search, is given in Dodds et al (2019), building on classical problem mapping techniques in Chancellor et al (2016, 2017). Secondly, it is impossible to map the problem Hamiltonian to qubits without specifying the solution outright. Hence, the search problem serves as a useful toy problem, especially in contexts where having analytic, computational, and physical implementations available for comparisons facilitates benchmarking and other testbed procedures.

The Sherrington–Kirkpatrick (SK) spin glass Hamiltonian H_SK (Sherrington and Kirkpatrick 1975) is defined on n spins as

$$\begin{eqnarray}&&{H}_{\mathrm{SK}}=-\displaystyle \frac{1}{2}\displaystyle \sum _{(j\ne k)=0}^{n-1}{J}_{{jk}}{S}_{j}{S}_{k}\end{eqnarray} \tag{ 15 }$$

where S_j are the classical spins (${S}_{j}\in \{-1,1\}$) and the couplings J_jk are drawn independently from the normal distribution ${ \mathcal N }(\mu ,{\sigma }_{\mathrm{SK}}^{2})$ with mean μ and variance ${\sigma }_{\mathrm{SK}}^{2}$. Finding the ground state of this Hamiltonian is NP-hard (Choi 2010), and uniformly hard, due to its finite-temperature phase transition (Sherrington and Kirkpatrick 1975).

It is computationally convenient to break the spin inversion symmetry by adding single-body field terms of the form ${\sum }_{j=0}^{n-1}{h}_{j}{S}_{k}$, where h_j are the field strength values. Like the couplings J_jk, the fields h_j are also drawn independently from ${ \mathcal N }(\mu ,{\sigma }_{\mathrm{SK}}^{2})$. When the fields strengths h_j are drawn from the same distribution as the coupling strengths J_jk, the hardness of finding the ground state follows directly from the hardness of the h_j = 0 case. The SK spin glass with such fields is mathematically equivalent to a zero field spin glass with one more spin which is 'fixed' in one orientation. This is not true in general for different distributions of field strength h_j. There are known examples in which fields can destroy spin glass behaviour (see, e.g. Young and Katzgraber 2004, Feng et al 2014). In particular, if the field strengths are much larger than the coupling strengths ($| {h}_{j}| \gg | {J}_{{jk}}|$ for all j, k), then the energy is minimized trivially when all the spins each minimize the energy with respect to their individual fields. While the distribution of field strengths could be used to tune the problem hardness, we do not use it in this way here, and only consider cases where the field and coupling strengths are drawn from the same distribution.

An astute reader will notice that if one effectively un-fixes the spin which corresponds to the fields (thus making all states two fold degenerate and converting the system to a double cover of the orignal system), these couplings will effectively be on average stronger by a factor of $\sqrt{2}$. As this increase in coupling strength does not scale with the number of spins, it is going to become less and less significant as the size of the system is scaled up the hardness will be preserved.

The mapping into the quantum Ising model is almost trivial: the classical spin variables S_j are simply mapped to Pauli-Z operators. Thus, the problem Hamiltonian ${\hat{H}}_{\mathrm{SK}}$ becomes

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{SK}}=-\displaystyle \frac{1}{2}\displaystyle \sum _{(j\ne k)=0}^{n-1}{J}_{{jk}}{\hat{Z}}_{j}{\hat{Z}}_{k}-\displaystyle \sum _{j=0}^{n-1}{h}_{j}{\hat{Z}}_{j},\end{eqnarray} \tag{ 16 }$$

The SK problem Hamiltonian differs from the search problem by having structure, produced by the ${\hat{Z}}_{j}{\hat{Z}}_{k}$ terms As a result, the covariances between the energies of two basis states depends on the Hamming-distance between them (Baldwin and Laumann 2018). Knowing the energy of one state gives some information about the energy of states that differ by a small number of bit-flips. This results in a distribution of the eigenenergies that is almost normal (as can be seen by plotting the distributions and numerically calculating moments), but which deviates from normal in the tails of the distribution.

To isolate the effect of the correlations in the SK problem, we compare it with the random energy model (REM) (Derrida 1980), in which the eigenenergies themselves are independently drawn from a normal distribution. The problem Hamiltonian ${\hat{H}}_{\mathrm{REM}}$ for REM is

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{REM}}=\displaystyle \sum _{j=0}^{N-1}{F}_{j}\left|j\right\rangle \left\langle j\right|,\end{eqnarray} \tag{ 17 }$$

with $\{\left|j\right\rangle \}{}_{j=0}^{N-1}$ the computational (Z) basis and the energies F_j drawn independently from the normal distribution ${ \mathcal N }(0,{\sigma }_{\mathrm{REM}}^{2})$.

REM has a similar energy level distribution to that of SK, apart from the tails. By definition it lacks the correlations: knowing the energy of one state gives no information about the energies of other states. Comparison between these two models highlights the effect of the pairwise structure in the SK model.

In contrast to the search problem, for SK and REM it is impossible to efficiently calculate the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ that maximizes the success probability. It is not even clear which measure of success probability should be maximized because, unlike the search problem, there will be no efficient way to find the optimal measurement time ${t}_{f}^{(\mathrm{opt})}$ for any choice of hopping rate γ. To bootstrap the investigation, we choose to define the optimal hopping-rate ${\gamma }_{\mathrm{opt}}^{(h)}$ with respect to one of the average probabilities defined in (9); in particular, we choose the hopping rate that maximizes the infinite-time average probability ${P}_{\infty }$ defined in (10). We make this choice because the infinite-time average probability ${P}_{\infty }$ is numerically convenient to calculate, and because it has been seen to be a relevant measure of probability in the search example, see figure 2(b). We will see in section 5.3 that the probability ${P}_{\infty }$ typically agrees well with probabilities averaged over shorter and more practical time windows.

Some plots of the infinite-time probability ${P}_{\infty }$ against hopping rate γ for typical 11-qubit examples of the SK and REM are shown in figure 3. Note that the maximal success probability varies by an order of magnitude between the two problem-types, with REM highest and SK lowest. While the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ is instance-dependent, these plots show that the dependence of infinite-time probability ${P}_{\infty }$ on hopping rate γ is typically characterised by broad, bumpy peaks for SK, and by narrow, well-defined peaks for REM. This implies that a precise value of the hopping rate γ is needed for REM, while there is some tolerance to non-optimal values of the hopping rate γ for SK for the sizes that we have studied.

Zoom In Zoom Out Reset image size

To investigate the success probability more systematically, we performed a brute-force numerical search to find the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ that maximizes the success probability ${P}_{\infty }$ for each spin glass instance from the data sets of 10k random instances for 5 ≤ n ≤ 11. This gives a baseline maximum average single run success probability for the quantum walk algorithm.

The optimal hopping rates ${\gamma }_{\mathrm{opt}}^{(h)}$ correspond to the best a quantum walk algorithm on the hypercube can possibly do in a single run. For practical algorithms, we need a heuristic method for choosing the hopping rate that can be calculated from the known parameters. For the quantum walk search algorithm, the optimal hopping rate balances the energy between the two components of the Hamiltonian, ${\hat{H}}_{P}$ and ${\hat{H}}_{G}$. Guided by this, we define the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ for SK and REM such that it balances these overall energy-scales on average. We match the energy-spread ${E}_{N-1}^{(h)}-{E}_{0}^{(h)}$ of the hypercube quantum walk Hamiltonian ${\hat{H}}_{h}$ with the average energy-spread $\langle {E}_{N-1}^{(P)}-{E}_{0}^{(P)}\rangle$ of the problem Hamiltonian ${\hat{H}}_{P}$. For the hypercube Hamiltonian ${\hat{H}}_{h}$ defined in (12), we have the energy spread ${E}_{N-1}^{(h)}-{E}_{0}^{(h)}=2n\gamma ;$ hence, we define the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ by

$$\begin{eqnarray}&&{\gamma }_{\mathrm{heur}}^{(h)}\equiv \displaystyle \frac{1}{2n}\left\langle {E}_{N-1}^{(P)}-{E}_{0}^{(P)}\right\rangle .\end{eqnarray} \tag{ 21 }$$

To demonstrate that this heuristic is sensible, we compare in figure 4 the distributions of optimal hopping rates ${\gamma }_{\mathrm{opt}}^{(h)}$ for SK (blue) and REM (red), as well as the heuristic hopping rates (black, dashed and dotted lines for SK and REM respectively) calculated according to (21), for the 11-qubit data set. For both SK and REM, the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ falls in the centre of the ${\gamma }_{\mathrm{opt}}^{(h)}$ distributions. Note that the SK distribution is much broader than for REM: not only are the individual peaks for ${\gamma }_{\mathrm{opt}}^{(h)}$ for SK much broader than for REM (figure 3), but the distribution of the maxima of those peaks is also much broader (figure 4). This may seem to be a problem for specifying a heuristic value for γ for SK from average energies, but as we will show, it is actually REM that fails for the heuristic γ, while SK works well.

Zoom In Zoom Out Reset image size

For a normal distribution of energy levels, the average problem energy spread can be estimated as

$$\begin{eqnarray}&&\left\langle {E}_{N-1}^{(P)}-{E}_{0}^{(P)}\right\rangle \simeq -({2}^{\tfrac{3}{2}}{\sigma }_{P}^{(\mathrm{energy})}){\mathrm{erf}}^{-1}\left(\displaystyle \frac{1}{N}-1\right),\end{eqnarray} \tag{ 22 }$$

where ${\sigma }_{P}^{(\mathrm{energy})}$ is the standard deviation of the energy eigenvalues of the problem Hamiltonian. For REM, the standard deviation ${\sigma }_{\mathrm{REM}}^{(\mathrm{energy})}$ is equal to the energy unit ${\omega }_{\mathrm{REM}}$ by definition (see section 3.2). For SK, the standard deviation ${\sigma }_{\mathrm{SK}}^{(\mathrm{energy})}$ can be shown to be equal to $\tfrac{{\omega }_{\mathrm{SK}}}{2}{[n(n+3)]}^{\tfrac{1}{2}}$. Equation (22) is accurate for REM (which has normally-distributed energy levels by definition) but, as already noted, the distribution of the eigenenergies in SK deviates from normal, especially in the tails. Numerically, we find that there is a multiplicative constant factor of approximately 0.887 that corrects the formula in (22) for SK for the effects of the non-normal tails. For the numerical analysis, we use the numerically calculated average energy-spread at each number of qubits n.

Figure 5(a) compares the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ and average optimal hopping rate $\langle {\gamma }_{\mathrm{opt}}^{(h)}\rangle$ at different numbers of qubits 5 ≤ n ≤ 11. The full width at half maximum (FWHM) has also been calculated for each instance, to estimate the tolerance ${\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}$ to deviations from the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ (illustrated in figure 3). The width of the shaded regions in figure 5(a) corresponds to the average tolerance range $\langle {\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}\rangle$ at each number n of qubits. While the heuristic hopping rate differs slightly from the the average optimal hopping rate for SK, the average tolerance range $\langle {\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}\rangle$ is much broader, and does not shrink with increasing number of qubits n. For REM, however, while we see close agreement on average, the tolerance range shrinks quickly with the number of qubits n as the peaks (as in figure 3, right) become narrower. This means that the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ is more likely to lie further than $2{\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}$ outside of the actual probability peak for each instance, even though it agrees well with the average optimal hopping rate $\langle {\gamma }_{\mathrm{opt}}^{(h)}\rangle$. Consequently, a quantum walk with the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ does not perform well for most REM instances.

Zoom In Zoom Out Reset image size

It is instructive to quantify this sensitivity to deviations from the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$. Figure 5(b) shows log–linear and log–log plots of the average fractional tolerance range $\langle {\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}/{\gamma }_{\mathrm{opt}}^{(h)}\rangle$ against number n of qubits for SK (blue circles), REM (red squares) and search (green triangles) on the hypercube. For SK, the fractional tolerance range $\langle {\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}/{\gamma }_{\mathrm{opt}}^{(h)}\rangle$ decreases as approximately 1/n, while for REM and search the decrease is approximately ${N}^{-0.5}$. This decrease is expected theoretically for search (Childs and Goldstone 2004). The fitted lines do not show exactly a square-root dependence (exponent of −0.5) due to the finite size effects for small numbers of qubits n ≤ 12.

Thus, we see that REM behaves like the search problem in a quantum walk setting. For a precisely optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$, the success probability is high, but this instance-dependent hopping rate is hard to predict, unlike for the analytically tractable quantum walk search algorithm. Without this precise hopping rate, quantum walks perform no better than guessing for the search problem and for REM. In contrast, quantum walks applied to SK give a better-than-guessing success probability ${P}_{\infty }\gt 1/N$ for the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ calculated according to (21).

With the conditions under which we can achieve a better-than-guessing success probability characterised for the three problem types, SK, REM, and search, we turn to the scaling of this success probability with problem size N.

Figure 6 shows how the single-time success probability P(t_f) varies with the measurement time t_f for two typical 11-qubit examples of SK and REM. In the REM case, the behaviour is similar to that shown in figure 2(a) for search: an oscillatory nature indicating the dominance of a two-level avoided-crossing feature, but with evidence of the population of other energy-levels that lead to finite-size effects in search. For REM, these finite-size effects are more pronounced, and are instance-dependent. The random nature of the REM problems means there is not such a clear cut off size, as there is for the search problem, above which finite size effects are negligible. In any case, based on search, we expect finite size effects to be significant at n = 11. For SK, the behaviour is quite different from search or REM. There is no indication of dominant oscillatory behaviour; instead, these plots show unpredictable, highly instance-dependent fluctuating dynamics for all the sizes we are using. This indicates that for SK, the behaviour is determined by the excitation of many energy levels.

Zoom In Zoom Out Reset image size

As with finding a suitable hopping-rate γ, both REM and SK differ from the search problem in that there is no practical way to find the optimal measurement time ${t}_{f}^{(\mathrm{opt})};$ a different approach must be taken instead. As already noted for the search problem, this can be handled by using the time averaged probabilities defined in (9). We first consider the infinite-time probability ${P}_{\infty }$, as defined in (10), since it is easy to calculate (see section 4). Figure 7 shows the average infinite-time success probability $\langle {P}_{\infty }\rangle$ against the number n of qubits for the two problems using both the optimal ${\gamma }_{\mathrm{opt}}^{(h)}$ and heuristic ${\gamma }_{\mathrm{heur}}^{(h)}$ hopping rates. For SK, this gives exponential decay with the number of qubits n in both cases: the average probability $\langle {P}_{\infty }\rangle$ changes with n according to

$$\begin{eqnarray}\langle {P}_{\infty }\rangle =\begin{array}{ll}\tilde{O}({N}^{-0.402\pm 0.001}) & {\rm{with}}{\gamma }_{\mathrm{opt}}^{(h)}\\ \tilde{O}({N}^{-0.417\pm 0.002}) & {\rm{with}}{\gamma }_{\mathrm{heur}}^{(h)}\end{array},\end{eqnarray} \tag{ 23 }$$

where $\tilde{O}$ may neglect factors logarithmic in its argument. That is, using the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ instead of the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ has only a minor impact on the average success probability $\langle {P}_{\infty }\rangle$.

Zoom In Zoom Out Reset image size

For REM, the behaviour is quite different. With the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ we see a success probability ${P}_{\infty }$ of constant order but with a pronounced dip. This behaviour is similar to that seen for the search problem, where the dip seen in figure 2(b) is a finite-size effect. This similarity is expected, given the similarity between the dynamical behaviour shown in figure 2(a) for search and in figure 6(b) for REM. With the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ for REM, we see a significantly reduced success probability ${P}_{\infty }$ compared to the optimal case. That is, the heuristic is performing poorly, despite the good agreement shown in figure 5(a).

The clear difference in behaviour between SK and REM can be explained by the different tolerances ${\rm{\Delta }}{\gamma }_{\mathrm{opt}}^{(h)}$ to deviations from the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ shown in figures 5(a) and (b). For SK, the tolerance range is broad enough for the heuristic to lie within it, while for REM the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ almost always misses this range entirely even though it is close to the average optimal hopping rate $\left\langle {\gamma }_{\mathrm{opt}}^{(h)}\right\rangle$.

We have thus numerically determined an average success probability scaling with problem size of $\sim \tilde{O}({N}^{-0.42})$ for a quantum walk finding SK spin glass ground states, using the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$. This is based on the infinite time-success probability ${P}_{\infty }$, i.e. uniform sampling from the distribution of all possible run times. We now investigate the time dependence in more detail: can we sample from a finite run time and still obtain the same speed up? Since P(0) = 1/N corresponds to random guessing, there must be a minimum time before which it is not effective to measure.

We define a mixing-time ${\tau }_{\mathrm{mix}}^{(\epsilon )}$ to be the latest time, t, for which the time averaged probabilities $\bar{P}(0,t)$ and $\bar{P}(0,2t)$ at the two times t and 2t differ by a fraction greater than the fluctuation parameter ,

$$\begin{eqnarray}&&{\tau }_{\mathrm{mix}}^{(\epsilon )}=\max \left\{t:\left|\displaystyle \frac{\bar{P}(0,t)-\bar{P}(0,2t)}{\bar{P}(0,t)}\right|\gt \epsilon \right\}.\end{eqnarray} \tag{ 24 }$$

This definition of ${\tau }_{\mathrm{mix}}^{(\epsilon )}$ is based on similar definitions found in prior work (Aharonov et al 2001), with modifications for computational convenience. We numerically estimated the mixing-time ${\tau }_{\mathrm{mix}}^{(0.05)}$ for each SK instance up to n = 11 qubits, using the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$ for each instance. We simulated the quantum walk computation dynamics for a successively-doubling duration until a time at which the condition is met was reached. The fluctuation parameter  = 0.05 corresponds to a deviation of 5%. To verify that the mixing-time ${\tau }_{\mathrm{mix}}^{(0.05)}$ correctly captures the relevant dynamical timescale, we also numerically estimated it for the search problem at each system size from n = 5 to n = 30 qubits. The search problem using continuous-time quantum walks can be mapped to the symmetric subspace, allowing larger sizes to be analysed. The mixing-time ${\tau }_{\mathrm{mix}}^{(0.05)}$ for search exhibits the expected exponential timescale: the solid green line of best fit in figure 8(a) has the expected scaling with problem size N of ${\tau }_{\mathrm{mix}}^{(0.05)}=\tilde{O}({N}^{1/2})$.

Zoom In Zoom Out Reset image size

For search, the scaling is dominated by the run time, the success probability is O(1). However, this behaviour only emerges clearly above n ∼ 20. Below this, the behaviour is influenced by the finite-size effects that arise due to population of higher energy levels. This means it is not useful to analyse the behaviour of the REM time scaling, finite size effects mask the scaling behaviour for computationally tractable sizes. However, unlike search and REM, the SK behaviour is influenced by higher energy levels at all sizes, through the frustration provided by the random couplings between the spins. Hence, we do not expect to see such finite-size effects in SK; the behaviour is already dominated by the frustration at small sizes. Figure 8(b) shows a log–log plot of the mixing-time (scaled by ω_SK) averaged over the ensemble $\langle {\tau }_{\mathrm{mix}}^{(0.05)}{\omega }_{\mathrm{SK}}\rangle$. The solid blue line of best fit has a logarithmic scaling with problem size N of

$$\begin{eqnarray}&&\langle {\tau }_{\mathrm{mix}}^{(0.05)}{\omega }_{\mathrm{SK}}\rangle =O({n}^{0.74\pm 0.03})\simeq O({[{\mathrm{log}}_{2}N]}^{0.75}).\end{eqnarray} \tag{ 25 }$$

Thus it contributes a logarithmic factor to the overall scaling. We emphasise that while this single-run timescale is polynomial in the number of spins n, the overall timescale is still exponential in n due to the exponential number of repeats required to achieve O(1) success probability.

To confirm the subsidiary nature of the time scaling for each SK run, we show in figure 9 a log-plot comparing, for the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$, the success probability ${P}_{\infty }$ in the infinite-time case (as in figure 7) and in the case of an early, logarithmically-scaling (with respect to N) measurement window $12.5{n}^{\tfrac{1}{2}}/{\omega }_{\mathrm{SK}}\leqslant t\leqslant 17.5{n}^{\tfrac{1}{2}}/{\omega }_{\mathrm{SK}}\equiv ({t}_{\mathrm{short}},{\rm{\Delta }}{t}_{\mathrm{short}})$. This n^0.5 scaling of the window is even shorter than the fitted scaling of n^0.75, although at these sizes the difference is not significant. This finite-time probability $\bar{P}({t}_{\mathrm{short}},{\rm{\Delta }}{t}_{\mathrm{short}})$ is similar to the infinite-time probability ${P}_{\infty }$: the solid purple line of best fit in figure 9 has a scaling with problem size N of

$$\begin{eqnarray}&&\left\langle \bar{P}({t}_{\mathrm{short}},{\rm{\Delta }}{t}_{\mathrm{short}})\right\rangle =\tilde{O}({N}^{-0.410\pm 0.002}).\end{eqnarray} \tag{ 26 }$$

This should be compared with (23), where the value of the exponent for the average infinite time success probability ${P}_{\infty }$ with the heuristic hopping rate ${\gamma }_{\mathrm{heur}}^{(h)}$ is given by −0.417 ± 0.002.

Zoom In Zoom Out Reset image size

As the dominant factor in the total runtime comes from the required number of repeats, and because the single-run timescale contributes only a logarithmic factor, these results constitute good numerical evidence for an average total runtime which scales with problem size N as $\sim \tilde{O}({N}^{0.41})$ for using quantum walks to find spin glass ground states, over the range of N in our data sets. This scaling is a better than the best possible (quadratic) speed up achievable for quantum walk search algorithms. Moreover, it comes without the requirement for exponential precision in setting the hopping rate that renders practical use of quantum walk searching difficult for large problems. We now present some insights into where the improvement over search comes from.

To investigate whether the energy correlations with Hamming distance in SK play a significant role in the computational process of finding the ground state with a quantum walk, we performed three additional sets of numerical tests.

Firstly, we used the same SK instances but performed the quantum walk using a complete graph Hamiltonian ${\hat{H}}_{K}$, defined in (11), instead of the hypercube graph Hamiltonian ${\hat{H}}_{h}$. This removes the correspondence of Hamming-distance between classical states with the distance between those states on the graph—for the complete graph, every state is one unit (edge) away from every other state. In terms of the Hamiltonian, the transverse Ising term is replaced by sums of products of up to n Pauli-X operators that flip up to n qubits at the same time, in all possible combinations. For each SK instance up to n = 11, we estimated the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(K)}$ for the complete graph, and then used it to calculate the infinite-time probability ${P}_{\infty }$.

Secondly, we constructed 'scrambled SK' instances, denoted sSK, by randomizing which state corresponds to which energy in the SK instances. In doing so, we arrive at Hamiltonians with identical energy spectra to the SK instances, but without the correlations between energy difference and Hamming distance on the hypercube graph. This approach has similarities with previous work (Farhi et al 2008, 2011, Hen 2014). For each sSK instance, we estimated the optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$, which is different from that used for the ordinary SK versions. This hopping rate was then used to calculate ${P}_{\infty }$.

Thirdly, we sorted the eigenenergies of each REM instance in increasing size and assigned them to the computational basis states in the order of a binary-reflected Gray code on their bitstrings, to arrive at a problem denoted REMGC. In doing so, we added some amount of Hamming-distance structure by ensuring that the closest energies are assigned to states that differ by only a single bit-flip. For each REMGC instance, we estimated an optimal hopping rate ${\gamma }_{\mathrm{opt}}^{(h)}$, which is different from that used for the ordinary REM problem. This was used to calculate the infinite-time probability ${P}_{\infty }$. While REMGC is not a hard problem as defined, it provides a useful example to compare with how the quantum walk finds the ground state of a SK spin glass.

These three variants provide separate tests of the influence of the graph structure (choice of quantum walk Hamiltonian) and problem structure (pairwise correlations in SK). Figure 10 shows how the infinite-time probability ${P}_{\infty }$ varies with the number of qubits n for these three variants, alongside SK and REM on a hypercube graph from figure 7. The variation of ${P}_{\infty }$ with the number of qubits n for the five variants is clearly split into two groups, behaviour like REM and search on the one hand, and behaviour like SK on the other. Removing the correlations from SK by scrambling the energies (sSK) results in behaviour like REM and search. Moreover, removing the correspondence between distance and Hamming weight by using the complete graph instead of the hypercube also changes the SK problem behaviour to be like REM and search. In the opposite direction, inserting pairwise correlations into REM via a Gray code (REMGC) results in problems that are much more like SK than like the REM problems on a hypercube graph.

Zoom In Zoom Out Reset image size

From this, we infer that the problem structure—in this case the pairwise correlations in SK—needs to be matched by a compatible driver Hamiltonian—in this case the hypercube/transverse Ising—to obtain better than quadratic scaling. This type of local structure in the solution space is exploited in many classical algorithms. For example, classical Monte Carlo optimizations that use a single bit flip update rule are naturally using this hypercube structure. Using a complete graph instead would correspond to flipping a random number of bits, which is equivalent to guessing at each step.

Continuous-time quantum walk time evolution is unitary, and there is no time dependence in the Hamiltonian that can lead to energy gain or loss by the system. Hence, it is important to consider how it can find a lower energy state than it starts in (with respect to ${\hat{H}}_{P}$) with any better-than-guessing probability. For the search problem, this happens through an analog of Rabi flopping (see figure 2), cycling between the initial and solution states. However, the dominant avoided level crossing structure is not present in the spin glasses to provide this mechanism.

We now show that there is a very generic mechanism (also described independently by Hastings 2019) that relies on starting in the ground state of the quantum walk part of Hamiltonian ${\hat{H}}_{G}$. Let $\langle \hat{O}{\rangle }_{\psi (t)}$ for operator $\hat{O}$ be defined by $\langle \psi (t)| \hat{O}| \psi (t)\rangle =\langle \hat{O}{\rangle }_{\psi (t)}$. Then, by linearity, and the definition of $\hat{H}(\gamma )$ in (8), the energy expectation at time t is

$$\begin{eqnarray}&&\langle \hat{H}(\gamma ){\rangle }_{\psi (t)}=\langle {\hat{H}}_{G}{\rangle }_{\psi (t)}+\langle {\hat{H}}_{P}{\rangle }_{\psi (t)}.\end{eqnarray} \tag{ 27 }$$

Due to the unitarity of the evolution under a time-independent Hamiltonian, this expectation energy will not change over time, giving

$$\begin{eqnarray}&&\langle \hat{H}(\gamma ){\rangle }_{\psi (t)}=\langle \hat{H}(\gamma ){\rangle }_{\psi (0)}.\end{eqnarray} \tag{ 28 }$$

which yields

$$\begin{eqnarray}&&\langle {\hat{H}}_{G}{\rangle }_{\psi (t)}-\langle {\hat{H}}_{G}{\rangle }_{\psi (0)}=\langle {\hat{H}}_{P}{\rangle }_{\psi (0)}-\langle {\hat{H}}_{P}{\rangle }_{\psi (t)}.\end{eqnarray} \tag{ 29 }$$

As $| \psi (0)\rangle$ is chosen to be the ground state of ${\hat{H}}_{G}$, the lhs must be non-negative. Furthermore, as $| \psi (0)\rangle$ is not an eigenstate of $\hat{H}(\gamma )$, some dynamics are guaranteed to occur and so the LHS must become positive at early times. Therefore, the rhs must also be non-negative always and positive at early times. Thus, taking any final time t_f, we get the inequality

$$\begin{eqnarray}&&\displaystyle \frac{1}{{t}_{f}}{\int }_{t=0}^{{t}_{f}}{\rm{d}}t\langle {\hat{H}}_{P}{\rangle }_{\psi (t)}\lt \langle {\hat{H}}_{P}{\rangle }_{\psi (0)}.\end{eqnarray} \tag{ 30 }$$

Equation (30) shows that performing time evolution under the computational quantum walk Hamiltonian from the initial state $\left|\psi (0)\right\rangle$ is guaranteed to lower the energy of the system with respect to ${\hat{H}}_{P}$ (the expectation value $\langle {\hat{H}}_{P}{\rangle }_{\psi (t)}$). This implies that the overlap with low energy eigenstates of ${\hat{H}}_{P}$ will increase, at least for short times. A measurement in the computational basis will thus be on average more likely than a random guess to produce a low energy state.

Starting in a low energy state is thus important for the success of the quantum walk algorithm (we have checked this numerically). It also implies that encoding prior information into the initial state will help, provided this is given in the form of a lower energy state than the uniform superposition state. This could be the final state from a previous run, for example, which will be explored further in Nita et al (2020). It is also necessary to bias the quantum walk Hamiltonian so that its ground state matches this biased initial state. Since this starting state is a known computational basis state, it is possible to do this biasing for suitably designed hardware.

For many optimization problem applications, it is helpful to find a low energy state, even if it is not actually the true ground state. From this point of view, that quantum walks necessarily lower the expectation energy with respect to the problem Hamiltonian is very appealing as a computational mechanism. This argument by itself does not provide a guaranteed scaling or quantum speed up, but it does explain how the quantum walk dynamics work in this setting, where there is no way to lose (or gain) energy. It is possible to generalise these arguments beyond time-independent Hamiltonians (Callison et al 2020a), to include monotonic functions A(t) and B(t) in (4).

To illustrate this energy redistribution mechanism, the plots in figure 11 show how the expectation value $\langle {\hat{H}}_{G}{\rangle }_{\psi (t)}$ of the quantum walk Hamiltonian (green solid-line) and the expectation value $\langle {\hat{H}}_{P}{\rangle }_{\psi (t)}$ of the problem Hamiltonian (red solid-line) vary during a quantum walk. We have included the instantaneous success probability P(t) (faint grey) to show that the timescale used is long enough for significant dynamics to take place. A typical 10-qubit SK example is shown in figure 11(a) and a typical 10-qubit REM example is shown in figure 11(b), both on the hypercube using their respective optimal hopping rates ${\gamma }_{\mathrm{opt}}^{(h)}$. Also shown is the ground state eigenvalue $\langle {\hat{H}}_{P}{\rangle }_{{E}_{0}^{\left(P\right)}}$ of the problem Hamiltonian (red, dashed–dotted line) and the ground state eigenvalue $\langle {\hat{H}}_{G}{\rangle }_{\psi (0)}$ of the quantum walk Hamiltonian (green, dashed line). In both SK and REM, the initial evolution takes the state away from the ${\hat{H}}_{G}$ ground state, raising the ${\hat{H}}_{G}$ expectation value, and thereby lowering the ${\hat{H}}_{P}$ expectation value to a point around which it fluctuates for the duration simulated. This clearly shows the energy redistribution mechanism at work, and the short time scale over which it appears.

Zoom In Zoom Out Reset image size