Continuous-time quantum walks for MAX-CUT are hot

By exploiting the link between time-independent Hamiltonians and thermalisation, heuristic predictions on the performance of continuous-time quantum walks for MAX-CUT are made. The resulting predictions depend on the number of triangles in the underlying MAX-CUT graph. We extend these results to the time-dependent setting with multi-stage quantum walks and Floquet systems. The approach followed here provides a novel way of understanding the role of unitary dynamics in tackling combinatorial optimisation problems with continuous-time quantum algorithms.

More recently CTQWs have been applied to a wider range of combinatorial optimisation problems, such as Sherrington-Kirkpatrick spin glasses and the random energy model [10].For both of these problems CTQWs were found to have a better scaling than Grover's algorithm.CTQWs have also been applied to MAX-2-SAT [11,12].The approach is relatively simple: a time-independent Hamiltonian is applied and the resulting distribution is sampled to try and find Robert J. Banks: robert.banks.20@ucl.ac.uk the solution or approximate solutions to the optimisation problem.Due to the exponentially growing Hilbert-space, previous work on combinatorial optimisation problems with CTQW has largely been limited to tens of qubits [10][11][12], with the exception of unstructured search.However, it has been shown CTQWs can always do better than random guessing [13].
In this paper we explore the performance of CTQWs on MAX-CUT, though we expect many of the arguments to hold for a wide range of combinatorial optimisation problems.In Sec. 3 we introduce CTQWs for MAX-CUT in more depth.By exploiting the effective low dimensionality of the evolution for short run-times we explore the performance of CTQWs in this limit, independent of problem size in Sec. 4.
The key take-home of this paper is that by assuming thermalisation we can analytically reason about the behaviour of CTQWs, using a phenomenologically motivated model discussed in Sec.5.2.CTQWs for optimisation being a relatively new and promising method for tackling combinatorial optimisation problems [10][11][12][13]9], this approach provides physical insight where exact diagonalisation is prohibitively difficult.
We emphasise all the discussion in this paper relates to CTQWs in the closed-system setting.Throughout this paper we adopt the convention h = 1 and k B = 1.The Pauli matrices are denoted by X, Y and Z.We also make use of the python packages QuTip [40,41] and NetworkX [42].

Introducing the Eigenstate Thermalisation Hypothesis
The aim of this section is to provide a very brief introduction to the ETH, focusing on the results pertinent to this paper.This section draws heavily from [16].For a more comprehensive review the reader is referred to the original work.This section is split into two parts: the first will discuss what is meant by a steady state, while the second part will introduce the consequences of the ETH.

The steady-state
Typically, classical random walks tend to steady states [43].Due to the oscillatory behaviour of the Schrödinger equation we generally do not expect this to be true for the state-vector of a CTQW.Instead of looking at the state-vector we might look at the expectation value of an observable O, with initial state |ψ i ⟩ and Hamiltonian H: where |E k ⟩ denotes an eigenstate of H.It follows that if the expectation of O is going to tend to a certain value, it will tend to its long timeaveraged value: From Eq. 1 it is clear that for the system to reach Eq. 2 there must be sufficient dephasing on sufficiently short time scales, such that the second sum in Eq. 1 becomes negligible compared to the first sum.In practice for any finite system we expect there to be fluctuations around the steady state value of O (i.e., around ⟨ Ō⟩).In short, a system is in its steady state if ⟨O(t)⟩ = ⟨ Ō⟩ for all t after a certain value, up to negligible fluctuations in ⟨O(t)⟩.For discussion about the time taken to reach equilibrium, please see [44].
Though Eq. 2 is useful, it requires complete knowledge of the spectrum of H, which quickly becomes intractable for large problem sizes.It is also perhaps not clear from this expression which observables correspond to steady states.The next section begins to address some of these questions.

The Eigenstate Thermalisation Hypothesis
The ETH is an attempt to explain how thermalisation occurs in a closed-system.The details can be found in [16].A consequence is that for a large class of observables the steady state value is equal to averaging over the microcanonical ensemble in the thermodynamic limit.The fluctuations from this value are exponentially suppressed with the degrees of freedom in the system [16].In certain cases a single eigenstate is sufficient to carry out the microcanonical averaging [15].
As mentioned prior, not all Hamiltonians exhibit ETH and not all observables thermalise [45].Typically, the Hamiltonian is assumed to be highly non-degenerate [15,14,24] and to provide a locally conserved quantity [18].The observables that exhibit thermalisation are typically considered to be local or few-body operators [15,24].For example, consider the highly non-local operator for a non-trivial system [18]: where E j ̸ = E k .This will clearly not approach a steady state.
In summary, the ETH tells us that for large systems the steady state, exhibited by local observables, can be well approximated by the microcanonical ensemble.Though there are known exceptions, this behaviour is thought to be very common [22,46].Note that it is possible for a system to approach a steady state that is not a thermal state [18].
For ease of computation, throughout this paper we make use of the canonical ensemble, which is equivalent to the microcanonical ensemble in the thermodynamic limit [24,28].That is to say, throughout this paper we assume the pure quantum state is locally well approximated by a thermal Gibbs state: with temperature 1/β.A CTQW is a time-independent Hamiltonian and in most settings unlikely to be integrable.We might therefore expect it to exhibit thermalisation for local observables.We provide numerical evidence for this in Sec.5.1.By exploiting the observation that the system is well approximated by a thermal state we make predictions on the performance of CTQWs.
3 Quantum walks for MAX-CUT

Set-up of the CTQW
In this section we introduce MAX-CUT for CTQWs.Given a graph G = (V, E), a cut separates the nodes into two disjoint sets.The weight of the cut is given by the number of edges which connect the two sets.The aim of MAX-CUT is to find the maximum cut.
This can be encoded as finding the ground state of an Ising Hamiltonian, In this paper we focus on graphs with fixed degrees and binomial graphs.To generate the binomial graphs each edge is selected with a certain probability, throughout this paper this probability is taken to be 2/3.Binomial graphs are also known as Erdős Rényi random graphs.
As the corresponding driver Hamiltonian we take the transverse-field Hamiltonian: where n = |V |.This corresponds to a quantum walk on the Boolean hypercube.
The CTQW Hamiltonian is given by: where γ > 0 is a parameter that needs to be set To assess the performance of Eq. 8 we look at is the ground state energy of H p , then the CTQW has found the ground-state.Random guessing would correspond to ⟨H p ⟩ = 0.The value of ⟨H p ⟩ can be related to the length of cut by: In short, the more negative ⟨H p ⟩ is, the better the performance of the CTQW.Energy is conserved during the evolution.Since the evolution starts in the ground-state of ⟨H d ⟩, ⟨H d ⟩ cannot decrease initially.Assuming non-trivial dynamics, this means that ⟨H p ⟩ must decrease to conserve energy.Hence, why Eq. 8 will do better than random guessing [13].
A second consequence of the conservation of energy is 0 is the ground-state energy of H d .This means that the evolution cannot completely occupy the groundstate of the CTQW Hamiltonian.This is in stark contrast to a number of other popular approaches to tackling optimisation problems with NISQ algorithms [47][48][49].
Finally, we note the Hamiltonian in Eq. 8 commutes with the spin-flip operator: Since, |+⟩ is an eigenstate of G with eigenvalue plus one, the evolution is restricted to the plus one eigenspace of G.This can be utilised in simulation to reduce the size of the Hilbert space by half.It also means that the CTQWs cannot reach a separable computational basis state, resulting in the ideal final state having non-zero entanglement.We also touch upon this in Appendix B in the context of the ETH.In the next section we discuss how to optimise γ to give the best possible performance.

Optimising the free parameter
The CTQW contains a free-parameter γ.Heuristically this controls the amount of dynamics present in the evolution.If γ is too large, we have approximate evolution under H p , leaving the measurement statistics unchanged.Correspondingly ⟨H p ⟩ ≈ 0. If γ is too small we have approximate evolution under H d , resulting in trivial dynamics and ⟨H p ⟩ ≈ 0. Clearly for a CTQW to be successful γ must be chosen somewhere between these two limits.
In this section we attempt to optimise the free parameter γ to give the best time-averaged value of ⟨H p ⟩.More explicitly the time-averaged value is given by: Numerically, to evaluate ⟨ Hp ⟩ we calculate: where |E k ⟩ is an eigenstate of H QW with corresponding eigenvalue E k .The intuition behind these equations can be found in Sec. 2.
In the rest of this section we find by brute force optimisation the optimal choice of γ for different graph choices.We also make comparisons to a reasonable heuristic choice of γ.To begin with we focus on a specific example of a 12-qubit binomial graph shown in Fig. 1, before looking at statistics gathered from multiple graph instances.
Starting with a specific choice of γ for the graph instance shown in Fig. 1, Fig. 2 shows a CTQW with γ = 1.The blue line in Fig. 2 shows the average value of H p from the Schrödinger equation, or more explicitly:  The dashed pink line shows the ground state energy of H p (i.e.E (p) 0 ).For the majority of the evolution shown ⟨H p (t)⟩ is fluctuating around the steady-state value, ⟨ Hp ⟩ (the dashed purple line).We can compare this to the ground state probability, shown in Fig. 3.The ground state probability shows significant oscillations, illustrating why we avoid optimising over this quantity.
In previous works, the free parameter in a CTQW has been heuristically chosen in an attempt to maximise dynamics.One approach to  do this has been to match the energy scales of the driver and problem Hamiltonian [10,11].For a MAX-CUT problem, this could be interpreted as Tr γ 2 H 2 p = Tr H 2 d such that: where κ 2 = |E| is the number of edges in the graph.
Focusing first on the performance of the binomial graph shown in Fig. 4, evaluating Eq. 14 gives γ heur = 0.5, corresponding to ⟨ Hp ⟩ min ≈ −5.88.This is a significant reduction in the performance of the CTQW.Fig. 5 shows the optimal γ for one hundred randomly generated problem instances for binomial graphs with problem sizes ranging between ten  and thirteen qubits.The optimal γ is seen to decrease with the problem size.The range of optimal values of γ is on the order of 0.1.The dashed red line represents the heuristic in Eq. 14 assuming the number of edges is n(n−1)/3.Clearly Eq. 14 is significantly underestimating the optimal γ.
As mentioned in the introduction of this section, it is expected that the optimal γ will balance H p and H d , hence why it is reasonable to assume γ ∝ n −1/2 .This is the same functional dependence on n as Eq.14.The dashed black line in Fig. 5 shows γ ∝ n −1/2 fitted to the available data.The curve is not inconsistent with the data.
Fig. 6 shows how ⟨ Hp ⟩ changes between γ opt  and γ heur for the same problem instances as Fig. 5.The difference in performance is on the order of one and appears to increasing with the problem size.There is clearly scope for improvement on the heuristic choice of γ.
We now turn to three-regular graphs, Fig. 7 shows γ opt for different problem sizes.For regular graphs the number of edges for a given problem size is fixed and scales with n.Hence, γ heur evaluates to be the same for all instances of a d-regular graph.This is the red dashed line in the figure .For three-regular graphs γ heur also appears to be generally underestimating γ opt but less so than in the binomial case.This is reflected in the change in performance, as shown in Fig. 8.
In this section we have focused on finding the optimal γ.The closer γ is to the optimal value, the better the approximation ratio is.We have also explored the use of a heuristic based around the mantra of maximising dynamics.The heuristic choice of γ has elucidated how the optimal choice of γ might scale with n.This is particularly pertinent to the discussion in Sec.5.2.In Sec.5.4 we show how the heuristic choice of γ (i.e.Eq. 14) can be recovered assuming thermalisation and how it might be improved upon.
An alternative, more involved, method for finding heuristic values of γ also based on the notion of maximising dynamics can be found in [13].
For the majority of this paper we are interested in the steady-state behaviour of the CTQW.However, there may be some advantage to running the CTQW on a very short time scale, with the aim of having a short time to solution.In the next section, we analytically investigate the very short time behaviour of CTQWs for MAX-CUT.

The very-short-time limit
A CTQW consists of a relatively simple, timeindependent Hamiltonian.Following [50,51] in this section we evaluate the torsion of the quantum evolution.The torsion is a local measure of how much the evolution deviates from a twodimensional space.By working in this sub-space we can make short-time predictions for ⟨H p ⟩.
The torsion depends only on expectations of powers of H QW , hence are constants over the duration of the evolution for a time-independent Hamiltonian.For the details of the calculation see Appendix A. The torsion of the evolution is given by: where κ 2 is the number of edges in the graph, κ 3 the numbers of triangles, and κ 4 the number of squares.
We can see that the torsion is reduced by the number of triangles in the graph.This suggests that the CTQW explores the solution space slower for graphs with a large number of triangles.The error from leaving a two dimensional space in a time ∆t can be approximated by [50]: So for times much less than T − 1 4 we can approximate the evolution of CTQWs by a two dimensional system.The details can be found in Appendix A. Under these assumptions, with γ > 0, the value of ⟨H p ⟩ can be calculated as: where As we can see, for very short times, ⟨H p ⟩ depends very little on the properties of the underlying MAX-CUT graph, depending primarily on the number of edges in the graph.The frequency term ω does depend on the number of triangles but for large problem sizes it is reasonable to expect ω 2 ≈ γ 2 κ 2 .Physically, it is reasonable, that at short times the CTQW only sees triangles and edges.Longer evolution would result in the approach seeing larger loops in the graph.
From a computational point of view this suggests that any useful short-time sampling must not be on a timescale much smaller than T − 1 4 since you would like the CTQW to see the whole graph.A rough estimate for a CTQW to see the whole structure of the graph might be lT − 1  4 , where l is the length of the largest loop in the problem graph.This is bounded from above by nT − 1 4 .In Fig. 9 we compare Eq. 17 to numerical simulation with a 12-qubit binomial graph with γ = 0.05.The dashed purple line shows the location of T − 1 4 , For short times there is good quantitative agreement between the numerical simulation (the blue line) and the two-level prediction (the pink line).At longer times there is reasonable qualitative agreement, with the twodimensional approximation capturing the oscillatory nature.
Examining a different 12-qubit binomial graph with γ = 1 gives Fig. 10.Again, we can see good agreement between the numerical simulation (the blue line) and the two-level prediction (the pink line) for short times.The dashed purple line corresponds to T − 1 4 with the dashed green line corresponding to nT − 1  4 .Here the two-level approximation provides poor qualitative insight into the CTQW outside of the short time limit, with the CTQW approaching an approximate steady state.Noticeably, by nT − 1 4 , the system has settled into the steady-state.
In this section we have explored the shorttime limit of CTQW for MAX-CUT.The argument has relied on a geometric property of the quantum evolution (i.e., the torsion).By integrating the Schrödinger equation in the reduced two-dimensional sub-space we have made a prediction on the performance of CTQWs at short times, without limiting the analysis to small problem sizes.For longer times, Eq. 17 breaks down.In general we are unable to integrate the Schrödinger equation for arbitrarily large problem sizes.For long times we turn to the Eigenstate Thermalisation Hypothesis (ETH).

Continuous-time Quantum walks in the Thermodynamic limit
As described in Sec. 2 a closed quantum system under evolution by a constant Hamiltonian can result in a state that is locally indistinguishable from a thermal state.In terms of CTQWs we are primarily interested in ⟨ Hp ⟩, a sum over local observables.Thus if a CTQW close to the optimal γ exhibits thermalisation, then ⟨ Hp ⟩ should correspond to that of a thermal state.In Sec.5.1 we provide numerical evidence that CTQWs do exhibit thermalisation even for small problem sizes.The aim of the rest of the section is to use this insight to estimate the optimal choice of γ, and ideally the associated value of ⟨ Hp ⟩.To do this, in Sec.5.2 a model for the density of states (DOS) is provided.The DOS is then used to predict the temperatures and entropy of entanglement associated with CTQWs in Sec.5.3.In Sec.5.4 we use the model to estimate ⟨ Hp ⟩ and to provide estimates for the optimal choice of γ.

Continuous-time quantum walks are well modelled by thermal states
Here we demonstrate, for a few examples, that CTQWs are well modelled by thermal states, taking ⟨ Hp ⟩ as a measure of success of the CTQW.
If we know the temperature associated with the CTQW, this provides us with a second route to find the same value of ⟨ Hp ⟩ by preparing a thermal state, thus side-steping the need for complete unitary dynamics.Since the Hamiltonian is constant during the evolution under a CTQW, energy is conserved in a closed-system.Hence the following must hold true for the thermal state [18]: where ρ β (H QW ) denotes a thermal state with Hamiltonian H QW and inverse temperature β.This equation fixes β, therefore there are no free parameters to fit.Fig. 11 shows ⟨H p ⟩ the Schrödinger evolution for a binomial graph and a three-regular graph and the corresponding prediction for ⟨H p ⟩ from the thermal state.For both problem instances it appears that ⟨H p ⟩ is fluctuating around a steady state value.Importantly, despite being far from the thermodynamic limit at only 12-qubits, the thermal state prediction is capturing the steady state behaviour well.Fig. 12 shows the solution to Eq. 19 for the optimal choice of γ for multiple instances of binomial and three-regular graphs.For both cases the inverse temperatures are quite small, with β < 1 for almost all instances.This means CTQWs correspond to Gibbs states with high temperatures.As mentioned in Sec. 2, we expect the ETH to hold true in the thermodynamic limit.Here we are working with of order ten qubits, hence we expect there to be finite-size effects meaning there will be some error between the Schrödinger equation and the thermal prediction.This is captured in Fig. 13.Generally, the thermal prediction, ⟨H p ⟩ β is overestimating the performance of the CTQW.As n is increased we would expect this error to decrease.However, even for these very small systems the error is relatively small, with no significant outliers in this data set.

Modelling the density of states
In the previous section we numerically demonstrated that even for relatively small systems, ⟨ Hp ⟩ is well predicted by a thermal state.Numerically solving Eq. 19 is however difficult, requiring finding the exponential of a matrix that increases exponentially with the problem size.In Sec.5.3 we will attempt to estimate the temperature associated with a CTQW.Before this is possible we need to model the density of states (DOS) for H QW .
To make progress we assume that the DOS is well modelled by a continuous distribution.We justify this assumption by arguing that the eigenstates of H QW are likely to become exponentially close as the system size grows.This follows from the largest eigenvalue of H QW being polynomial in n, and there being exponentially many states.Assuming no (or little) degeneracy, the difference between energy levels must shrink exponentially with the problem size.We also assume that the DOS is well modelled by a uni-modal distribution.This assumption will break down if γ is too large or too small.However, as seen in Sec.
3.2 we expect useful values of γ to correspond to somewhere between these two limits.
Intuitively, for large systems, we would expect the energy of a randomly chosen eigenstate to be close to the average energy.The average energy, µ, of H QW is given by: Now consider the variance, σ 2 : Therefore, as n goes to infinity, σ/2 n tends to zero.So we expect a tightly peaked distribution for large n.The details for calculating moments of the DOS can be found in Appendix C. The DOS for thermalising systems has previously been modelled with a Gaussian distribution [52].Since the CTQW exhibits thermalisation we adopt this approach.Fig. 14 shows a Gaussian fit to the DOS for a 12-qubit three-regular graph (Fig. 14a) and a binomial graph (Fig. 14b) with γ optimised for each problem.Visually this appears  to be an acceptable approximation for the regular graph, less so for the binomial graph.Given the simplicity of this approximation we utilise it throughout the paper to make simple analytic predictions.
To verify the Gaussianity of the DOS, we look at two typical tests of normality: the skewness, s, of the distribution and the excess kurtosis k.The skewness is the ratio of the third moment of the distribution to the cube of the standard deviation.The kurtosis is the ratio of the fourth moment of the distribution to the variance squared.The excess kurtosis is the kurtosis minus the expected kurtosis of a Gaussian distribution, which is equal to 3.
Starting with the skewness: Further details can be found in Appendix C. Note that triangle-free graphs have no skewness and that for positive γ, the skewness is always nonnegative.
For a d-regular graph, κ 2 = dn/2 and κ 3 is bounded by O d 2 n .So the skewness, following the heuristic in Eq. 14 (provided γ is held constant) scales approximately as n −1/2 .Therefore it will tend to zero as the problem size is increased.
Conversely, for binomial graphs, provided γ scales proportional to 1/ √ n, then the skewness will not scale with n.Hence, the Gaussian approximation will not hold as well for a binomial graph.
Examining now the excess kurtosis gives: Under the same assumptions above for the scaling of the skewness, the excess kurtosis will tend to zero as n tends to infinity for regular graphs.For binomial graphs the excess kurtosis is unlikely to vanish.
For large regular graphs the Gaussian approximation looks to hold well.The same cannot be said to be true for binomial graphs.This is a consequence of regular graphs looking locally tree-like in the infinite sized limit [53].Hence it may be necessary to consider other models for the DOS for binomial graphs and for regular graphs away from the infinite-size limit.
Here we propose using an exponentially modified Gaussian (EMG) distribution [54,55] to model the DOS for a CTQW.The EMG is a skewed unimodal distribution.The skew can only be positive, as is the case with the DOS for a CTQW.In the correct limit it can recover a Gaussian.It has the following probability density function: where erfc is the complementary error function.
The fitting parameters m, ν and λ can be ascer-  tained from n, σ and s.For the details the reader is referred to Appendix D.2.The Gaussian approximation only includes information about the problem-size and therefore will make the same prediction for a large number of graphs.For instance, all regular graphs with the same degree have the same DOS under the Gaussian approximation.By incorporating the skewness into the model we are incorporating more information about the graph structure, namely the number of triangles in the problem.Visually, as shown in Fig. 15, we can see that the EMG distribution (the dashed red line) provides a better fit than the Gaussian distribution for the binomial graph.The EMG distribution also still models the DOS for the regular graph well.During the rest of this paper we use the DOS to make analytic predictions about the behaviour of a CTQW in a closed-system setting.

Predictions from the density of states
Having provided a model for the DOS for H QW , in this section we make two predictions in relation to CTQWs, namely the temperature and the entropy of entanglement.

Estimating the temperature
Finding the temperature of the CTQW requires finding the solution to Eq. 19.To find approximate solutions we use the DOS models to evaluate the partition function, Z. Eq. 19 then be-comes [56]: Assuming a Gaussian DOS then the estimated inverse temperature is given by: This is perhaps what one would estimate as the temperature on dimensional grounds alone.The full details can be found in Appendix D.1.
Taking the EMG to model the density of states gives: where σ 2 = n + γ 2 κ 2 and ∆ = γ(3κ 3 ) 1/3 .The details can be found in Appendix D.2.Fig. 16 shows the error in temperature between β est and β (the numerical solution to Eq. 19).The approximate solutions are consistently underestimating the inverse temperature.The error for binomial graphs is substantial, with typical errors being between 25% and 30% (Fig. 16a).For three-regular graphs, this reduces to somewhere between 15% and 25% (Fig. 16b).Using the EMG distribution (i.e.Eq. 25) provides little improvement on the estimate for the threeregular graphs (Fig. 17b).This is unsurprising given how close to Gaussian the DOS is for these problems.For binomial graphs the improvement is more substantial with with typical errors being between 10% and 15%.As n is increased, we expect that the DOS will be better modelled by a continuous DOS, hence we would expect the error in β to be improved.
By assigning the quantum evolution a temperature, we have mapped the challenge of understanding a dynamical problem to a static problem.So far in this section we have shown how to reasonably estimate the associated temperature.If the temperature is too high, then the associated thermal state can be efficiently classically approximated.Results by Crosson et al. [57] suggest that for values of β ≤ 0.1 a classical computer could simulate the associated thermal state efficiently for a three-regular graph, suggesting that for a CTQW to provide an advantage, it must operate outside this regime.

Estimating the entropy
Often considered of high interest for a quantum algorithm is the entanglement in the system [58].This can be quantified by the entanglement entropy.Here we calculate the entanglement entropy by tracing out half of the qubits, which are randomly selected, and calculating the von Neumann entropy of the reduced density operator.
We compare the entanglement entropy to the thermodynamic entropy derived from the modelled DOS, given by: Since entropy is an extensive quantity we take half of the above and compare it to the entanglement entropy.Although the ETH cannot be straightforwardly applied here as the entropy is not a local observable, the two have been shown to be linked in previous works [59][60][61].
For the Gaussian distribution Eq. 26 reduces to: and for the EMG: Fig. 18 shows the entanglement entropy averaged over a hundred 14-qubit binomial graphs as a function of time.As is clear from the figure the entanglement entropy is approximately constant.The dashed purple line shows S norm and the dashed pink line S EM .Both are overestimating the entanglement entropy.However, they appear to be reasonable estimates with S EM providing the better estimate.As mentioned in Sec.
3, the ideal final state will be a low-entanglement state.CTQWs lack any mechanism to dissipate entanglement, as shown by Fig. 18 5.4 Estimating ⟨ Hp ⟩ In Sec.5.3.1 we saw how the temperatures associated with a CTQW could be found.By assuming the system is well modelled by a Gibbs distribution, it follows that holds approximately true.The derivation behind this equation can be found in Appendix E. For the Gaussian approximation this gives: which is optimised by The details can be found in Appendix D.1.Note that this is the same result as the heuristic of maximising the dynamics (i.e.Eq. 14).The same approach can be taken for the EMG model.The corresponding result is too cumbersome to be written out in full.Fig. 19 compares ⟨ Hp ⟩ to the predictions from the Gaussian and EMG predictions for a single 12-qubit binomial graph.Though both overestimate the performance of the CTQW, both provide reasonable estimates on the optimal γ.The EMG clearly provides a better model for ⟨ Hp ⟩ than the Gaussian model.Indeed it is remarkably close to the direct numerical calculation of ⟨ Hp ⟩.Note that although the thermal model appears to provide good estimates for ⟨ Hp ⟩ for large and small γ, in this regime ⟨H p (t)⟩ will not necessarily display steady-state behaviour.For example in Fig. 9, γ is very small and there are significant oscillations present.However, this regime of γ, as discussed in Sec. 3, is unlikely to be of practical interest to CTQWs.For the performance of γ opt est the reader is referred to Sec. 3.For the EMG distribution there are two important quantities to examine: • How much does the performance of ⟨ Hp ⟩ change between the optimal γ and the γ predicted by the EMG DOS?
• What is the difference between the prediction of ⟨ Hp ⟩ and the true value?
The first of these is addressed for small problem sizes in Fig. 20.The performance is improved over γ heur from simply balancing the drive and problem Hamiltonians, especially for binomial graphs.Hence, we believe we have found a good heuristic method for finding γ for CTQWs applied to MAX-CUT.Fig. 21 shows the error in the performance between the true value of ⟨ Hp ⟩ and the predicted value from the EMG approximation, for the γ predicted to give the best possible performance from the EMG DOS.The error is relatively small for these small problem sizes.Given the tractability of Eq. 29, and evidence of small errors, it is possible to estimate the performance for larger problem sizes.This is shown in Fig. 22.
By treating closed-system CTQWs as thermalising systems we have shown they thermalise.From here we have been able to provide useful heuristic choices for optimising γ and to make practical predictions for ⟨ Hp ⟩.For further discussion of the performance of CTQWs see Appendix F. We have also provided an alternative approach to understanding CTQW that might be applicable to other combinatorial optimisation problems.
6 Multi-stage quantum walks

The set-up
So far we have focused on time-independent Hamiltonians.In this section we introduce timedependence into the CTQW, by implementing multi-stage quantum walks (MSQWs).In this case γ is increased monotonically in a step-wise fashion.With each step of γ the system is allowed to reach equilibrium.For instance an l-stage quantum walk would have the following schedule: This schedule can be seen as being motivated by approaches like Quantum Annealing [47].For each step, while γ is held constant energy is conserved.Therefore we can modify Eq. 19 to be: QW is the Hamiltonian during the k th stage of the quantum walk (i.e. with γ k ) and |ψ (t)⟩ is the state-vector associated with the MSQW.Fig. 23 shows a five-stage CTQW for a 12-qubit binomial graph.Each stage increases γ by 0.5, starting with 0.5.Each stage has a duration of twenty units of time.With each stage ⟨H p ⟩ (the solid blue line) improves and quickly tends to an approximate steady-state.Solving Eq. 33 to find the temperature and associated performance gives the purple line in the figure.For this instance, we have good qualitative and reasonable quantitive agreement between the thermal prediction and the Schrödinger equation.During the first stage of the MSQW, the thermal prediction is quite far from the prediction of the Schrödinger equation, suggesting the state is far from thermal.This is to be expected as γ is small and the driver Hamiltonian dominates, breaking the assumption that the DOS is continuous.Despite this, as γ is increased the system thermalises.Perhaps most interestingly, as γ is further increased the state remains thermal despite the problem Hamiltonian becoming more dominant.The inverse temperatures associated with this MSQW are β = 1.21, 0.72, 0.53, 0.41, 0.32.This corresponds to heating the system, suggesting that despite the final stage providing the best performance, it might be the distribution easiest to simulate classically [57].Given the thermal behaviour, it might be reasonable to assume that the techniques developed earlier in this paper for CTQW can be applied to MSQW.

Making predictions
In Sec. 5 the tools for making predictions about CTQWs were developed.Here we apply these tools to the MSQW.Here we only use the EMG distribution to model the DOS, since this also captures the Gaussian model.
To make predictions Eq. 33 needs to be approximated to find β.Since finding |ψ(t)⟩ is likely to be numerically intractable for large systems we The difference between the stationary state calculated from the Schrödinger equation and a thermal state for five-stage CTQW on 100 12-qubit binomial graphs."EMG DOS" refers to the prediction with the energy fixed by Eq. 34.The numerical approach gives the prediction from using the numerically determined energy, numerically fixing the temperature and calculating ⟨H p ⟩.
In Fig. 24 we consider the final value of ⟨ Hp ⟩ for a hundred five-stage CTQW on 12-qubit binomial graphs compared to the analytical prediction (i.e.Eq. 34) and the numerical prediction assuming thermalisation.The schedule is the same as the schedule used in Fig. 23.Consider first the numerical prediction.Assuming thermalisation, the correct final steady state for all the problem instances is predicted within an error of 0.5.In contrast the EMG DOS model only achieves this error bound for approximately 50% of the instances.This approach suffers from cumulative errors in ⟨H p ⟩ at each stage of the MSQW (e.g.Eq. 34), which can result in large errors in the prediction.Here we have numeric evidence for MSQW exhibiting thermalisation.The EMG DOS model struggles to capture the performance for all instances, as well as the time-independent case.This is due to the difficulty in determining the energy of the system.To improve this approach a better model for the DOS or prediction of ⟨H p ⟩ needs to be developed, especially for the small γ case.

Implications for gate-based implementation
One approach for realising a CTQW might be to simulate it on a gate-based quantum computer.The simplest approach to achieve this is to Trotterise the evolution [62], such that the Hamiltonian alternates between H H p .The resulting state-vector is prepared by the quantum computer: where t = N s τ and N s is is an integer.Here N s acts a proxy for circuit depth.Setting τ presents a trade-off between minimising gate-depth with a large τ , and accuracy of the simulation that requires τ to be small.Alternatively, Eq. 35 could be viewed as a QAOA [63] schedule with each step having the same set of angles.In the previous sections the Hamiltonians have been held constant for extended intervals in time.This has provided a local conserved quantity, resulting in a finitetemperature thermal state.The Trotterised system is not time-independent and energy is not conserved.In this section we investigate Eq. 35.
Since the Hamiltonian is periodic in time, it is an example of a Floquet system (a very brief introduction can be found in Appendix G).
To analyse the long time behaviour of the Floquet system we evaluate: where |ϕ(t) α ⟩ are the Floquet modes, ϵ α the corresponding quasi-energies, and c α = ⟨ϕ(0) α |+⟩ [64].This is illustrated in Fig. 25 for a 12 qubit binomial graph with τ = 0.1 and γ = 1.The solid blue line shows the Schrödinger equation.Within a few intervals of τ the system has already thermalised to ⟨ Hp ⟩ F = −7.74.This is shown by the dashed purple line (i.e. the time averaged value from Eq. 36).The performance of the Floquet system is marginally worse than the value of ⟨ Hp ⟩ = −7.79 for the corresponding CTQW.Indeed Fig. 26 shows the difference between ⟨ Hp ⟩ F and the corresponding value for a CTQW for one-hundred binomial qubit instances.
For the Floquet instances τ = 0.2.For each problem instance the CTQW provided the better performance.
Depending on the size of τ , we might expect two things to occur for ⟨ Hp ⟩ F [65].For small τ the system does not have time to thermalise, so |ψ F (t)⟩ mimics the behaviour of a CTQW.For large τ , the system thermalises at each step.Since energy is not conserved, the system tends to the infinite temperature Gibbs state, resulting in ⟨ Hp ⟩ F = 0.This is shown by the blue line in Fig. 27 for a single 12-qubit binomial graph as τ is increased.

The thermal model
In this section we provide an approximate timeindependent model to capture the stroboscopic behaviour of the Floquet system assuming the Trotter-step is small.Previous works have focused on using a Magnus expansion [66][67][68], for ease of computation here we take an alternative approach.
The stroboscopic behaviour is captured by the single period unitary: from here it is a trivial exercise to recover the associated Hamiltonian: where U p (t) = e −iHpγt/2 .Assuming τ is small, making the approximation gives a time-independent Hamiltonian.We can also neglect the constant scaling of the Hamiltonian that is not going to change the long time behaviour.Hence we are left with the following time-independent Hamiltonian: From here we can follow the steps outlined in Sec. 5 for the CTQW.Note that Tr Hk so we can use the moments for the DOS derived for the CTQW.Therefore any change in performance must come from the energy of the initial state: where deg(k) is the degree of the k th node.Notably, for binomial graphs (or any graph where deg(k) increases with n) as n goes to infinity, ⟨+| HF |+⟩ goes to zero.This corresponds to an infinite temperature thermal state.Consequently, in the thermodynamic limit, under the above assumptions, any attempt to approximate a CTQW with a Floquet system will fail.In contrast, regular graphs will have more stability.The numerical prediction for thermalising under HF can be seen in Fig. 27 (dashed purple line) for a single 12-qubit binomial graph.The numerical prediction agrees well with the true numerical value of ⟨ Hp ⟩ F (solid blue line) up until the crossover to an infinite temperature Gibbs state.
As established in the earlier section on CTQWs, they correspond to Gibbs states.Therefore, we use the same model as for the CTQW with a modified energy constraint (i.e.Eq. 42).This is shown by the dashed pink line in Fig. 25, for a single 12-qubit binomial graph to good agreement.The results for one hundred 12-qubit binomial graphs with τ = 0.2 and γ = 1 can be seen in Fig. 28.As shown in the figure the median error between the numerically calculated value of ⟨ Hp ⟩ F and the analytic value, assuming an EMG DOS and thermalisation, is close to zero.
Perhaps most usefully, this thermal model provides some intuition on how mitigate Trotter errors for this system.Essentially, we wish to start in the effective ground-state of the Floquet system.Eq. 40, suggests this might be achieved by changing the initial state from |+⟩ to U † p (τ /2) |+⟩.This is illustrated for a single 12-qubit binomial graph in Fig. 29 with γ = 1.The dashed purple line shows the steady-state value for a CTQW.The solid blue shows the steady state of the Floquet system with the initial state being the |+⟩ state.The pink line shows the same Floquet system with the modified initial state U † p (τ /2) |+⟩.Indeed, for the short time steps, the modified initial state better matches the CTQW.Indeed it even provides a slight advantage for a small interval in τ .

Conclusion
Throughout this paper we have attempted to understand the performance of CTQWs outside of what is classically simulable for non-integrable models.For the short-time limit, we were able to exploit the effective low dimensionality to demonstrate that the performance can be characterised by underlying graph properties, such as the number of edges in the graph.
For the steady-state we conjectured that the system thermalises.This means despite the statevector consisting of 2 n complex amplitudes, for a given γ, the value of ⟨H p ⟩ will depend on one real number, the energy.This provides some insight into the computational mechanisms involved in CTQWs.Assuming thermalisation, classical statistical physics provides an alternative route to understanding CTQWs for MAX-CUT on a broad range of graphs within a closed-system setting.By associating the unitary dynamics with a temperature, it provides an alternative route to achieve the same value of ⟨ Hp ⟩ , either through dissipative dynamics or classical simulation.Here we have introduced an EMG distribution to account for some of the frustration in the model to make analytic predictions.By exploiting this model we were able to find reasonable estimates for the optimal choice of γ, that utilised properties of the underlying graph.Importantly, we were able to make predictions far away from what it is easy to directly and completely simulate.We have also demonstrated some preliminary results for MSQWs.
Further to this we were able to extend these results to Floquet systems as an approach to implementing CTQWs.For all the problems considered the Floquet system performed worse than the corresponding CTQW with the same γ.Our results predict that in the thermodynamic limit, MAX-CUT on binomial graphs will thermalise to an infinite temperature Gibbs state, independent of the step size.This is consistent with known work on thermalisation in Floquet systems [65].The Floquet systems considered in this paper can be thought of as a QAOA [63] ansatz with restricted angle choices.For QAOA we wish to avoid the system thermalising to an infinite temperature Gibbs state, hence it would be prudent to avoid a Floquet ansatz.Further to this it is reasonable to assume large step sizes might be detrimental for high-depth QAOA circuits.It is has been observed that a QAOA circuit with randomly chosen step sizes also thermalises to an infinite temperature state [58], suggesting that thermalisation is not solely the consequence of the discrete time symmetry of the Floquet system.How far the results of this paper might extend to other NISQ algorithms such as QAOA and Quantum Annealing is an ongoing area of research.
A Details for the short time limit

A.1 Calculating the torsion
Here we provide the details for the short time limit for CTQW.Firstly, we calculate the torsion, and the related quantity curvature.The curvature is a local measure of how much the quantum evolution deviates from a geodesic.These are conserved quantities valid for the whole evolution.We then show how this leads to short time predictions on ⟨H p ⟩ for these problems in Sec.A.2.
From Laba et al. [50] the curvature of the wave-function, for a time-independent Hamiltonian is given by: where ∆H = H − ⟨H⟩.The torsion is given by: Therefore, calculating the torsion and curvature reduces to evaluating ⟨∆H j ⟩, for j = 2, 3, 4. Since these are all conserved quantities, we can take the expectation with respect to the initial state, |+⟩.The rest of this section details the rather tedious steps required to find these expectation values.Those willing to take the expression for T at face-value are invited to skip to Sec.A.2.The first step is to find ⟨H⟩: since each Z i Z j in H p will flip qubits i and j to the minus state.Consequently, only H d contributes to this expectation value.Consider now: where κ 2 = |E| is the number of edges in the graph.It was pointed out by Gnatenko et al. [51], that ⟨H 2 p ⟩ is equal to the number of edges in the graph.To see why this is the case, we can write out this term explicitly to give: where the sum is over the ordered pairs (i k , j k ), for k = 1, 2. The only possible way a term in the sum can be non-zero is if i 1 = i 2 and j 1 = j 2 .Therefore the number of non-zero terms corresponds to the number of edges in the graph.The rest of the terms in ⟨∆H 2 ⟩ are relatively trivial to calculate.Extending the above logic, Gnatenko et al. [51] demonstrated that: where κ 3 is the number of triangles in the graph, and where κ 4 is the number of squares in the graph.
Using the above results we can calculate ⟨∆H 3 ⟩ and ⟨∆H 4 ⟩.Starting with ⟨∆H 3 ⟩: Calculating ⟨H 3 ⟩ gives: where we have neglected to write out terms that are trivially zero and made use of Eq. 45.To find The X k in each term is not going to change which terms are non-zero.It will only introduce a minus sign if k is i 2 or j 2 .Hence ⟨H p H d H p ⟩ is equal to: Alternatively, we can argue as follows: in the absence of the X k term, Eq.48 would give −nκ 2 .Accounting for the edge cases where k = i 2 or k = j 2 gives: Combining all of the above we are left with: To find C and T , it remains to find ⟨∆H 4 ⟩.
Evaluating ⟨H 4 ⟩ gives: Again the presence of X in the terms of the sum is not going to change which terms evaluate to be non-zero.Fixing the edge (i 2 , j 2 ), if either one of the X terms coincides with this edge and the other X term does not, it contributes −1 to the sum; there are 4(n−2) ways of this happening.If both X terms do not coincide with this edge then the term contributes 1; there are (n − 2) Finally, we can assemble all the expectation values to find The resulting expression for the curvature is: and for the torsion: Note that the torsion is zero for a two-qubit system, with one edge, as expected from the spin-flip symmetry in the problem.The torsion is also zero for a ring of three qubits.

A.2 Short-time estimate for the performance
In the previous section we calculated the torsion.The torsion measures how much the wave-function deviates from a two-dimensional space in a time ∆t.The error from deviating from this space is given by [51]: with ε 2D = 0 corresponding to remaining in the two-dimensional subspace and ε 2D = 1 corresponding to departing the subspace.If ε 2D ≪ 1, then the two-dimensional approximation holds well.That is for times t ≪ T −1/4 , we can apply a two-level approximation.At t = 0 the state of the system is |ψ (t = 0)⟩ = |+⟩.Integrating the system forward some infinitesimal time δt, gives |ψ (t = δt)⟩ = (I − iHδt) |+⟩.
Applying the Gram-Schmidt procedure to the two vectors, |ψ (t = 0)⟩ and |ψ (t = δt)⟩, gives the orthonormal basis for the subspace: where the expectation is with respect to the |+⟩ state.Writing out H in this subspace is: which evaluates to neglecting the uninteresting term proportional to the identity.The problem Hamiltonian, within this space is: evaluating this gives: where sgn (γ) denotes the sign of γ.
Having now explicitly evaluated the relevant operators in this subspace, it remains to calculate the expectation of the problem Hamiltonian, which is: where B Spin-flip symmetry and the ETH In Sec. 3 we pointed out MAX-CUT has a spin-flip symmetry.To recap the notation, given a CTQW for MAX-CUT with Hamiltonian H, then [H, G] = 0, where: Typically, the system starts in the |+⟩ state so the evolution is restricted to the plus one eigenspace of G for the entire evolution.Despite this restriction, in Sec. 5, we used the 'complete' Hamiltonian as opposed to the Hamiltonian projected onto the correct symmetry sector to apply the ETH.This was numerically verified in Sec.5.1.We justify this choice on the following grounds [18]: 1. G is a global operator and corresponds to a global symmetry.

We are interested in calculating the expectation values of local interactions only.
It is well established within the ETH literature that local symmetries can prevent thermalisation [18].In contrast the spin-flip symmetry is a global symmetry.In general, there are a large number of globally conserved quantities for evolution under a constant Hamiltonian.If H belongs to a Hilbert Space H with dimension dim H, then there are at least dim H globally conserved quantities [18].To see this, let the projector onto each eigenstate of H be denoted by P k with k = 1, . . ., dim H, then [18]: Hence ⟨P k ⟩ are globally conserved quantities.Despite the large number of conserved quantities, we generally do not expect these global symmetries to have much impact on the dynamics.Though the above argument suggests we should be careful when considering the role of globally conserved quantities, it fails to consider the specific case of the spin-flip symmetry which we know restricts dynamics to the positive eigenspace of G for a CTQW.Here we emphasise again that we are only interested in local quantities.A locally measurable consequence of the spin-flip symmetry is that: where we have used the initial state being an eigenstate of G and [U, G] = 0, consequently ⟨Z i (t)⟩ = 0.
Since ⟨Z i (t)⟩ is a local observable, this should agree with the result if we were to replace the unitary evolution with a Gibbs distribution: As [H, G] = 0 this is also zero.Thus the thermal-state recovers the correct expectation value for any local observable that flips sign under conjugation by G (i.e., GO L G = −O L ).That is not to say that the spin-flip symmetry has no role in calculating expectation values.Consider the case where the initial state is a linear superposition of two states, where If O L also commutes with G then (as is the case with H p ) it follows that: However, we have conjectured each (and numerically investigated) that each of the above matrix elements can be represented by replacing each unitary evolution with a thermal state.That is to say, therefore the system would need to be assigned two temperatures for each symmetry sector.
Throughout this paper we have tacitly assumed we are in the paramagnetic phase, justified by the associated high temperatures.Work on spontaneous symmetry breaking and the ETH can be found in [30,69].

C Calculating the moments of the density of states
To better understand the density of states (DOS) it is possible to calculate moments of the distribution.We can use these as fitting parameters to the model for the DOS.We are interested in calculating the moments of the distribution produced by the eigenenergies of the CTQW Hamiltonian for the graph G = (V, E): where: where n = |V | is the number of nodes in the graph.Note that the Hamiltonian is constructed of Pauli matrices, which are traceless [62].
Denoting the eigenenergies of H QW as E k , then the mean of the distribution of eigenenergies is given by: and since the Hamiltonian is traceless, this evaluates to zero.Repeating this approach for the variance gives: where κ 2 is the number of edges in the graph.Between the penultimate and final line, we have used that only terms that are equal to the identity will contribute to the trace.The same approach can be used to find the third moment: where κ 3 is the number of triangles in the graph.Further details can be seen in [51] or Appendix A. The fourth moment also follows from the above logic: From here it is straightforward to calculate the skewness and excess kurtosis.In theory, this process of moments could be continued to higher orders, incorporating loops of greater lengths.However, as shown above, even by the fourth order this becomes cumbersome.
then the exponentially modified Gaussian distribution is given by: Through application of the convolution theorem [71] it is straightforward to write down the moment generating function (and hence partition function) from the moment generating functions of the exponential and Gaussian distributions.The resulting partition function in terms of the fitting parameters m, ν 2 and λ is: To find the fitting parameters we fit the moments of h(x) to the moments associated with H QW .
Utilising the moment generating function to find the mean (µ), variance (σ 2 ) and skew (s) of h(x) gives: The above equations can be inverted to find the fitting parameters in terms of the mean, variance and skew associated with H QW .These properties can be extracted from H QW by calculating Tr H k QW for k = 1, 2, 3.The method can be found in Appendix C. From here the same steps as Appendix D.1 can be followed to derive predictions for the temperature, entropy, ⟨H p ⟩ and optimal choices of γ E Derivation of Eq. 29.
Using Eq. 29 removes the need for matrix exponentiation, hence here we discuss in greater detail the assumptions behind this equation.We assume that the system is well modelled by a Gibbs state and consider the function gives: Acting on both sides with e −βH QW and taking the trace gives:   ).The results can be seen in Fig. 30.For the small problem sizes considered the performance seems largely independent of system size.This perhaps is unsurprising for a high temperature Gibbs state with a local Hamiltonian [73,74].We might expect different parts of the system to be weakly correlated, hence the CTQW is optimising locally and does not scale with problem size.
In practice one is perhaps less interested in how the average of the energy distribution of ⟨H p ⟩ compares with the absolute minimum, particularly if considering CTQWs as an exact solver.In such a case one might be more interested in the ground-state probability or time-to-solution.This requires some notion of run-time not discussed in this paper and will be the subject of future investigations.

G A brief introduction to Floquet systems
In this appendix we provide a very brief introduction to Floquet systems, with the aim of providing sufficient background to understand Sec. 7.This introduction draws from [64].
In this paper we have used Floquet system to refer to a closed quantum system with a time periodic Hamiltonian.A closed system will evolve under the Schrödinger equation: We are interested in Hamiltonians that are time-periodic, such that: where τ is the period of the Hamiltonian.The Floquet theorem then tells us that there exist solutions to the Schrödinger equation of the form where ϵ α are termed quasi-energies and |ϕ α (t)⟩ the Floquet modes.The Floquet modes are time periodic, |ϕ α (t)⟩ = |ϕ α (t + τ )⟩.The quasi-energies are time-independent but are only uniquely defined up to multiples of 2π/τ .This is reminiscent of Bloch's theorem [75] for spatially periodic crystals.
Hence |ϕ α (t)⟩ is an eigenvector of a single period unitary with eigenvalue e −iϵατ .Therefore, by studying the single period unitary we can construct the wave function at all times, including the long-time average in Eq. 36.

Figure 2 : 1 .
Figure 2: A time domain plot of ⟨H p ⟩ (the blue line) for a randomly generated 12-qubit instance shown in Fig. 1.The dashed pink line shows the minimum energy of H p .The dashed purple line is the time-averaged value, according to Eq. 12.

Figure 3 :
Figure 3: A time domain plot of the ground-state probability for the randomly generated 12-qubit instance shown in Fig. 1.

Figure 4 :
Figure 4: The time-averaged value ⟨H p ⟩ as γ is changed for the 12 qubit graph shown in Fig. 1.The dashed purple line shows the location of γ opt and ⟨ Hp ⟩ min .

Figure 5 :
Figure 5: The optimal value of γ for one hundred problem instances for each problem size for binomial graphs.The red dashed line reflects the heuristic from Eq. 14, explicitly γ = 3/(n − 1).The dashed black line corresponds to fitting the curve γ opt = an −1/2 to the medians of the data, yielding a ≈ 2.93.

Figure 6 :
Figure 6: The difference in ⟨ Hp ⟩ for γ opt and γ heur for the same binomial graphs as Fig. 5.

Figure 7 :
Figure 7: The optimal value of γ for three-regular graphs.In ascending problem size, there are: 16, 45 problem instances.The red dashed line is the heuristic from Eq. 14, explicitly γ = 2/3.

Figure 8 :
Figure 8: The difference in ⟨ Hp ⟩ for γ opt and γ heur for the same three-regular graphs as Fig. 7.

Figure 9 :
Figure 9: The performance of a CTQW on a 12 qubit binomial graph with γ = 0.05.The blue line shows the result of direct integration of the Schrödinger equation.The pink line shows the result of the two-dimensional approximation (i.e, Eq. 17) The dashed purple line shows the location of T − 1 4 .

Figure 10 :
Figure 10: The performance of a CTQW on a 12 qubit binomial graph with γ = 1.The blue line shows the result of direct integration of the Schrödinger equation.The pink line shows the result of the two-dimensional approximation (i.e, Eq. 17) The dashed purple line shows the location of T − 1 4 .The dashed green line shows the location of nT − 1 4 , with n = 12.

Figure 11 :
Figure 11: The performance of a CTQW with γ optimised to give the best value of ⟨ Hp ⟩.The dashed purple line shows the thermal state prediction.The temperature is fixed using Eq.19.

Figure 12 :
Figure 12: The inverse temperature associated with CTQWs (i.e. the solution to Eq. 19).For each problem instance γ is equal to γ opt , shown in Fig. 5 and Fig. 7.

Figure 14 :
Figure 14: Histogram for the DOS for two graphs with optimised γ.The energies have been binned into 100 bins.The DOS has been normalised such that the total density is equal to one.The dashed black line shows the fitted Gaussian distribution.

Figure 15 :
Figure 15: Histogram for the DOS for two graphs with γ optimised.The energies have been binned into 100 bins.The DOS has been normalised such that the total density is equal to one.The dashed black line shows the fitted Gaussian distribution.The dashed red line shows the EMG distribution.

Figure 16 :
Figure16: The error in the predicted inverse temperature assuming a Gaussian DOS (i.e.Eq. 24).For each problem instance γ is equal to γ opt , shown in Fig.5and Fig.7.

Figure 17 :
Figure 17: The error in the predicted inverse temperature assuming an EMG DOS (i.e.Eq. 25).For each problem instance γ is equal to γ opt , shown in Fig. 5 and Fig. 7.

Figure 18 :
Figure 18: The entanglement entropy averaged over a hundred 14-qubit binomial instances (solid blue line).The dashed purple line shows S norm and the dashed pink line S EM The shaded regions show a single standard deviation.

Figure 19 :
Figure 19: The solid blue line shows ⟨ Hp ⟩ for a 12-qubit binomial graph.The dashed purple (pink) line shows the prediction from Eq. 29 assuming a Gaussian (EMG) DOS.

Figure 20 :
Figure 20: The difference in ⟨ Hp ⟩ between the optimal choice of γ and the prediction from the EMG DOS (γ EM ).

Figure 21 :
Figure 21: The difference between the prediction of ⟨H p ⟩ with an EMG DOS using Eq.29 (⟨H p ⟩ EM ) and the true value.For each instance the value of γ has been chosen such that ⟨H p ⟩ EM is minimised.Each bar shows one hundred binomial graph instances.

Figure 22 :
Figure 22: Estimated scaling of the CTQW for binomial graphs assuming thermalisation and an EMG DOS.Each data point shows the median performance of one hundred binomial graph instances.

Figure 23 :
Figure 23: A five-stage CTQW on a 12-qubit binomial graph.The purple line is the numerical prediction assuming that the MSQW is well modelled by a CTQW.The pink line shows the prediction from Eq. 29 and modelling the DOS as an EMG distribution.

Figure 24 :
Figure 24:  The difference between the stationary state calculated from the Schrödinger equation and a thermal state for five-stage CTQW on 100 12-qubit binomial graphs."EMG DOS" refers to the prediction with the energy fixed by Eq. 34.The numerical approach gives the prediction from using the numerically determined energy, numerically fixing the temperature and calculating ⟨H p ⟩.

Figure 25 :
Figure 25: A Floquet system simulating a CTQW for a 12-qubit binomial graph, with γ = 1 and τ = 0.1.The dashed purple line shows the time averaged value (i.e.Eq. 36).The dashed pink line shows the prediction from Eq. 29, assuming an EMG DOS.

Figure 26 :
Figure 26: The reduced performance resulting from implementing the CTQW as a Floquet system for onehundred 12-qubit binomial graphs.For the Floquet systems τ = 0.2.For all the problem instances γ = 1.

Figure 27 :
Figure 27: The steady state value for a single 12-qubit binomial graph with γ = 1 as τ is increased.The blue line shows the prediction from Eq. 36.The dashed purple line shows the prediction using the time-independent Hamiltonian shown in Eq. 40.

Figure 28 :
Figure 28: The difference between the true steady state and the steady state predicted assuming the system is well modelled by a time-independent Hamiltonian with a EMG distribution.The figure shows the results for one hundred 12-qubit binomial graph instances with γ = 1 and τ = 0.2.

Figure 29 :
Figure 29: Improved performance from changing the initial state.The dashed purple line shows the steady-state value for a CTQW.The solid blue line shows the steady state of the Floquet system with initial state |+⟩.The solid pink line shows the steady state of the Floquet system with initial state U † p (τ /2) |+⟩.

Tr 1 )
∂ γ e −βH QW = − QW , H p e −βH QW .(96)Evaluating the terms in the sum for which k > 1:TrH (k−1) QW , H p e −βH QW = j ⟨E j | H (k−1) QW , H p e −βH QW |E j ⟩ = j ⟨E j | H QW , H (k−2) QW , H p |E j ⟩ e −βE j = j E j ⟨E j | H (k−2) QW , H p − H (k−2) QW , H p |E j ⟩ e −βE j = 0,where |E j ⟩ denotes the eigenvectors of H QW .Therefore, once we assume the state is well modelled by a Gibbs distribution it follows that:Tr ∂ γ e −βH QW = −β Tr H p e −βH QW .(97)Commuting through the trace with the partial derivative and dividing both sides by the partition function gives Eq. 29.

Figure 30 :
Figure 30: The performance (i.e., ⟨ Hp ⟩) of the CTQW for the optimal γ predicted by the EMG DOS compared to the ground state energy of the problem Hamiltonian 1.The system is prepared in the ground-state of H d , here denoted by |+⟩, and evolved under H QW .The resulting state, after a time t, is 2 d H 2 p + H d H p H d H p + H d H 2 p H d + H p H 2 d H p + H p H d H p H d + H 2 Again, we have neglected to write out the terms that evaluate to zero.Most of the terms in this expansion we have already evaluated, or are simple to evaluate.We are left with ⟨H p H 2 d H p ⟩, ⟨H p H d H 2 p ⟩, ⟨H 2 p H d H p ⟩, and ⟨H 4 p ⟩, the last of which can be evaluated with Eq. 46.Turning now to ⟨H p H 2 d H p ⟩: 2ways of this happening.Finally if both edges coincide, this edge contributes 1; there are four ways of this happening.The result is:⟨H p H 2 d H p ⟩ = (n − 4) 2 κ 2 .Thus, it remains to evaluate ⟨H p H d H 2 p ⟩ and ⟨H 2 p H d H p ⟩.Note that: ⟨H p H d H 2 p ⟩ * = ⟨H 2 p H d H p ⟩.To evaluate ⟨H p H d H 2 p ⟩, consider ⟨H 3 p ⟩.As before, the introduction of H d is not going to change which terms in the sum are non-zero, if we were to write out a similar sum for this term, as in Eq. 48.Therefore, ⟨H p H d H 2 p ⟩ is going to depend on the number of triangles in the graph.The introduction of H d is going to result in sign flips to some terms.The result is Tr H 4 d + 4γ 2 H 2 d H 2 p + 2γ 2 H d H p H d H p + γ 4 H 4 p = n 2 + 2(n 2 − n) + 4γ 2 nκ 2 + 2γ 2 (n − 4)κ 2 + γ 4 (κ 2 + 3κ 2 (κ 2 − 1) + 24κ 4 ) .