Quantum Sampling Algorithms, Phase Transitions, and Computational Complexity

Drawing independent samples from a probability distribution is an important computational problem with applications in Monte Carlo algorithms, machine learning, and statistical physics. The problem can in principle be solved on a quantum computer by preparing a quantum state that encodes the entire probability distribution followed by a projective measurement. We investigate the complexity of adiabatically preparing such quantum states for the Gibbs distributions of various classical models including the Ising chain, hard-sphere models on different graphs, and a model encoding the unstructured search problem. By constructing a parent Hamiltonian, whose ground state is the desired quantum state, we relate the asymptotic scaling of the state preparation time to the nature of transitions between distinct quantum phases. These insights enable us to identify adiabatic paths that achieve a quantum speedup over classical Markov chain algorithms. In addition, we show that parent Hamiltonians for the problem of sampling from independent sets on certain graphs can be naturally realized with neutral atoms interacting via highly excited Rydberg states.


I. INTRODUCTION
A sampling problem is the task of drawing samples from an implicitly defined probability distribution, which may, for example, be the Gibbs distribution of a classical system at a fixed temperature. This particular class of sampling problems encompasses many practically significant and challenging problems with applications in areas as varied as statistical physics [1], optimization [2], and machine learning [3]. Moreover, sampling problems involving Gibbs distributions represent a natural setting to explore connections between phase transitions and computational complexity. It has been shown that Markov chain algorithms with local updates can be used to efficiently sample, in a time polynomial in the system size, from the Gibbs distribution of lattice models if and only if spatial correlations are short ranged [4,5]. This implies that it is easy to sample from many Gibbs distributions at high temperature, while no efficient, general-purpose algorithm is known below the ordering transition temperature of the underlying physical system.
The development of programmable quantum systems raises the question of what connections exist between quantum complexity and phase transitions. In analogy to the statement about classical Gibbs distributions, it has been shown that a quantum computer can efficiently generate samples from quantum Gibbs distributions with short-ranged correlations [6,7]. Quantum mechanics also offers an alternative motivation to consider sampling problems. Any quantum state |ψ together with an orthonormal basis {|s } encodes a sampling problem: According to the Born rule, a projective measurement in the {|s } basis yields the outcome s with probability p(s) = | s|ψ | 2 . By defining p(s) in terms a quantum gate sequence [8] or an optical network [9], it is possible to construct sampling problems that can be efficiently solved on a quantum computer but which are hard to solve classically. These efforts have led to impressive experimental demonstrations [10,11], showing that current quantum devices may outperform classical devices for specifically tailored sampling problems.
In our recent work, reference [12], we introduced a family of quantum algorithms that provide unbiased samples by preparing a state that encodes the entire Gibbs distribution. We showed that this approach leads to a speedup over a classical Markov chain algorithm for several examples. In this Article, we explore in detail a novel connection between classical sampling problems and quantum phase transitions that arises in this context. Concretely, we construct a so-called parent Hamiltonian H q , which has the defining property that the state |ψ encoding the probability distribution of interest is its ground state. The construction starts from a Markov chain that converges to the desired probability distribution. The gap between the ground state and the first excited state of the parent Hamiltonian can be related to the mixing time of the Markov chain, thereby establishing a connection between quantum phases and the classical complexity of the sampling problem. We also investigate the quantum complexity of the problem by exploring the time required to adiabatically prepare the state |ψ . We focus on sampling from classical Gibbs distributions at inverse temperature β, which naturally define a continuous family of sampling problems as well as a continuous family of states |ψ(β) and Hamiltonians H q (β). Analyzing in detail several examples, we find for each that adiabatic state preparation along the one-parameter family H q (β), which may be viewed as a quantum analogue of simulated annealing [13], exhibits the same time complexity as sampling by means of the Markov chain used to construct the parent Hamiltonian. By contrast, a quantum speedup can be achieved along paths in a suitably extended parameter space for all of the examples. In each case, we explain the magnitude of the speedup in terms of the distinct quantum phases occupying the extended parameter space and the nature of the phase transitions between them.
Our results complement existing quantum algorithms for sampling problems [13][14][15][16][17][18][19][20][21][22][23][24][25][26] by offering a physical picture of quantum speedup. Moreover, our approach is suitable for implementation on near-term quantum devices. We show in particular that the parent Hamiltonian for sampling from independent sets on certain graphs has a natural implementation using highly excited Rydberg states of neutral atoms that is compatible with recently demonstrated programmable atom arrays [27,28].
This Article is organized as follows. In section II, we show how the parent Hamiltonian H q (β) can be constructed from a Markov chain. In section III, we discuss general aspects of adiabatic state preparation. This formalism is applied to the first example, the Ising chain, in section IV. Other examples, namely sampling from independent sets and the unstructured search problem, follow in section V and section VI. Section V also includes a scheme to experimentally realize the parent Hamiltonians that arise for sampling from independent sets. We provide a summary and outlook in section VII.

II. PARENT HAMILTONIANS
We consider the problem of sampling from the Gibbs distribution of a classical Hamiltonian H c at inverse temperature β. Labeling the microstates of the system by s, the probability distribution is given by p(s) = e −βHc(s) /Z, where Z = s e −βHc(s) is the partition function. The entire Gibbs distribution can be encoded in the quantum state which we will refer to as a Gibbs state. For concreteness, we focus on systems composed of n classical spins such that s may be represented by a string of n bits. It is worth noting that the Gibbs state can be represented by a projected pair entangled state (PEPS) with bond dimension D = 2 under the additional assumption that the Hamiltonian involves only single-and two-body terms [29]. Preparing a quantum system in the Gibbs state allows one to obtain an independent sample from the Gibbs distribution associated with H c by simply performing a projective measurement in the {|s } basis. Hence, sampling from a classical Gibbs distribution can be reduced to preparing the corresponding Gibbs state. We can establish a connection between the complexity of state preparation and quantum phases by introducing the notion of parent Hamiltonians. A parent Hamiltonian of a state |ψ is any Hamiltonian for which |ψ is a ground state. Starting from another Hamiltonian whose ground state can be easily initialized, the Gibbs state can be efficiently prepared adiabatically if the initial Hamiltonian can be smoothly deformed into the parent Hamiltonian without encountering a small energy gap. By contrast, adiabatic state preparation is slow if the Hamiltonian exhibits a small gap (vanishing in the thermodynamic limit) at any point along the adiabatic path. This occurs in particular when the physical system undergoes a quantum phase transition. The parent Hamiltonian of a Gibbs state is not unique. If a PEPS representation of the state is known, there exists a general method of constructing a parent Hamiltonian [29,30]. For our purposes, it is more convenient to pursue a different construction, where the parent Hamiltonian is obtained from a Markov chain whose stationary distribution is given by the Gibbs distribution [29].
Markov chains constitute a powerful tool to solve sampling problems. To sample from the Gibbs distribution p(s) = e −βHc /Z, we introduce the generator matrix M , whose matrix elements M (s, s ) specify the transition probability from microstate s to s . A probability distribution q t (s) at time t evolves into q t+1 (s) = s q t (s )M (s , s) after one time step of the Markov chain. We impose detailed balance on the Markov chain, which can be expressed as e −βHc(s ) M (s , s) = e −βHc(s) M (s, s ). Detailed balance ensures that p(s) is a stationary distribution of the Markov chain, that is, p(s) = s p(s )M (s , s). If, in addition, the Markov chain is fully mixing and aperiodic, the Perron-Frobenius theorem guarantees that p(s) is the unique stationary distribution [31]. Hence, running the Markov chain for a sufficient number of steps yields a sample from the desired Gibbs distribution. The relevant number of steps is known as the mixing time, which can be related to the spectral properties of the generator matrix. Since M is a stochastic matrix, its largest eigenvalue is λ 1 = 1 with p(s) being the corresponding left eigenvector. Under the assumptions that render the stationary distribution unique, the largest eigenvalue is nondegenerate and the eigenvalue with the second largest magnitude satisfies |λ 2 | < 1. It can be shown that the mixing time satisfies the lower bound τ m ≥ |λ 2 |/(1 − |λ 2 |) [32,33]. This inequality shows that a Markov chain mixes slowly if the spectral gap between the largest and second largest eigenvalue is small.
The generator M can be turned into a parent Hamiltonian of the Gibbs state in Eq. (1) by a similarity transformation following reference [29] and related earlier work [34][35][36]. As a consequence of detailed balance, quantum Hamiltonian. The Gibbs state |ψ(β) is an eigenstate of H q (β) with eigenenergy zero, which follows from the fact that M is a right stochastic matrix. We can further show that |ψ(β) is the unique ground state of H q (β) by noting that for every eigenvalue n(1 − λ) of H q (β), there exists a corresponding eigenvalue λ of M . Since the largest eigenvalue λ 1 = 1 of M is nondegenerate, the ground state of H q (β) is also nondegenerate and has eigenenergy zero. The factor of n in Eq. (2) ensures that the spectrum of the parent Hamiltonian is extensive. The rescaled spectrum of the parent Hamiltonian as compared to the generator matrix M of the Markov chain reflects the natural parallelization in a quantum system, where each qubit represents a physical resource. For a fair comparison between the mixing time of the Markov chain and the adiabatic state preparation time studied below, we divide the mixing time by n, defining t m = τ m /n. It then follows from the correspondence of the spectra of M and H q (β) that t m ≥ 1/∆(β) − 1/n, where ∆(β) is the gap between the ground state and first excited state of the parent Hamiltonian.
The family of Hamiltonians H q (β) can be viewed as generalized Rokhsar-Kivelson Hamiltonians, which are stoquastic and frustration free [35][36][37][38]. The quantum phases of H q (β) and the phase transitions separating them not only allow us to understand the performance of the Markov chain in physical terms, but they also explain why adiabatic state preparation along certain paths in parameter space offer a quantum speedup over the Markov chain. The achievable speedup depends on the nature of the quantum phase transitions encountered along the adiabatic path and its origin can be traced to quantum phenomena such as ballistic transport and tunneling.

III. ADIABATIC STATE PREPARATION
While there exist general approaches to running Markov chains on a quantum computer with a guaranteed quantum speedup [16,17,23,24], these methods require deep circuits that are beyond the reach of current quantum devices. For this reason, we explore the prospect of preparing Gibbs states adiabatically. This alternative avenue has the potential of being realized on near-term devices if the adiabatic evolution can be naturally implemented with little overhead and low noise. Since the adiabatic state preparation time strongly depends on the local rate of change of the Hamiltonian, we now describe a heuristic scheme used to determine the rate of change along a fixed path in parameter space.
The scheme is based on the adiabaticity condition, which states that the population of the instantaneous eigenstates |n due to nonadiabatic transitions from the eigenstate |m can be approximately bounded by where E m and E n are the respective instantaneous eigenenergies [39,40]. We require that the transition probability from the ground state into any excited state be small along the adiabatic path. Denoting the ground state by |0 and setting the ground state energy E 0 = 0, we require the sum of the right-hand side of Eq. (3) over all excited states to satisfy for some small constant ε. For a fixed path in parameter space, Eq. (4) determines the rate of change of the Hamiltonian parameters. It is convenient to introduce a dimensionless parameter s, which increases monotonically from s = 0 at the beginning of the path to s = 1 at the final point. The Hamiltonian is parametrized by a set of parameters {λ µ }. An adiabatic path is fixed by the functional dependence λ µ (s), while the rate of change of the parameters is captured by the differential equation which follows immediately from Eq. (4), where and |∂ µ 0 is a shorthand for ∂|0 /∂λ µ . The total evolution time can be obtained by integrating Eq. (5), yielding where The above adiabatic schedule exhibits a number of important features. Equation (6) accounts not only for the gap between the ground state and the first excited state, but also the gap to all other excited states. In this way, we ensure that the adiabatic condition can be satisfied at second-order quantum phase transitions, where an extensive number of states is energetically close to the ground state. In addition, Eq. (6) takes into account the matrix elements that induce nonadiabatic transitions to appropriately weigh the significance of each excited state. The total evolution time t tot can be tuned by adjusting the constant ε. The constant of proportionality l only depends on the path and not its parametrization. For this reason, we refer to l as the adiabatic path length and to g µν as the adiabatic metric.
To analyze the performance of adiabatic state preparation, we introduce the fidelity where |ψ(β) is the desired Gibbs state and |φ(t tot ) is the final state after evolving under the time-dependent Hamiltonian for a given adiabatic path with the rate of change of the parameters as described above. The fidelity bounds the total variation distance d = ||p − q|| between the desired Gibbs distribution and the prepared distribution q(s) = | s|φ(t tot ) | 2 by d ≤ √ 1 − F [41]. The adiabatic theorem guarantees that F → 1 in the limit t tot → ∞. Here, we are interested in the shortest evolution time t a such that the fidelity exceeds some threshold, taken to be 1 − 10 −3 throughout. We will see below that t a corresponds to a small value of ε that is approximately independent of the system size. Hence, the adiabatic path length can serve as a proxy for the total time required to prepare a Gibbs state with high fidelity. One could thus in principle identify good adiabatic paths by minimizing the adiabatic path length, that is, by finding geodesics under the metric g µν . However, for the present purposes, it suffices to compare a small number of paths that uncover distinct asymptotic dependencies of the computation time on the system size.

A. Parent Hamiltonian and quantum phases
We now illustrate the formalism outlined above with several concrete examples, starting with the ferromagnetic Ising model in one dimension. To construct the parent Hamiltonian, we choose Glauber dynamics as the Markov chain, which updates the configurations according to the following prescription [42]: Pick a spin at random and draw its new orientation from the Gibbs distribution with all other spins fixed. Starting from configuration s, the probability of flipping spin i in the Ising chain is thus given by By promoting the values of the spins s i to operators σ z i , we can concisely write the generator of the Markov chain as Equation (2) then gives where 4h(β) = 1+1/ cosh(2β), 2J 1 = tanh(2β), 4J 2 (β) = 1 − 1/ cosh(2β) (see [43,44] for early derivations of this result). While Eq. (12) defines a one-parameter family of Hamiltonians dependent on β, it is natural to consider the quantum phase diagram of H q in an extended parameter space with arbitrary values of h, J 1 , and J 2 . The Hamiltonian in Eq. (12) is exactly solvable using a Jordan-Wigner transformation to map it onto free fermions (see appendix A 1 for details) [45,46]. One finds that the energy of an excitation with momentum k is proportional to At a quantum phase transition, the energy gap between the ground state and the first excited state vanishes, which implies that there must exist at least one mode with zero energy, i.e., E k = 0 for some k. This is indeed the case at k = 0 when h + J 1 + J 2 = 0 (line 1 in Fig. 1), k = π when h − J 1 + J 2 = 0 (line 2 in Fig. 1), or at k = π − cos −1 J1 2J2 when h − J 2 = 0 and |J 1 | < 2|h| (line 3 in Fig. 1). To identify the distinct phases, we consider line cuts through the parameter space. Along J 2 = 0, the Hamiltonian in Eq. (12) reduces to the transverse-field Ising model, which exhibits a well-know phase transition from a paramagnetic phase to a ferromagnetic phase at |J 1 | = |h|. Similarly, the line J 1 = 0 has been studied previously, allowing us to identify the remaining third phase as a cluster-state-like phase separated from the paramagnetic phase by a symmetry-protected topological phase transition [47,48].
All three quantum phases are gapped and the gap only vanishes at the phase transitions. The dispersion relation given by Eq. (13) is linear at low energies for all phase transitions with the exception of the tricritical points at J 1 /h = ±2, J 2 /h = 1. At the tricritical points, the low-lying excitations from two distinct phase transitions merge into a single gapless mode with a quadratic dispersion relation. This behavior is captured by the dynamical critical exponent z [49], which takes the value z = 2 at the tricritical points and z = 1 everywhere else. For a finite-sized system of n spins, the gap closes as ∼ n −z since the spacing of momenta is inversely proportional to n. We note that at the transition between the paramagnetic and cluster-state-like phases, the gap displays an oscillatory behavior as a function of system size and may even vanish exactly. Nevertheless, the envelope follows the expected ∼ n −1 scaling.
The above insights into the phase diagram and the nature of the phase transitions have immediate implications on both the mixing time of Markov chain and the adiabatic state preparation time. The parent Hamiltonian H q (β) (red curve in Fig. 1) occupies the paramagnetic phase and is therefore gapped for all finite β. Hence, the Markov chain is expected to mix rapidly in this regime as confirmed by the rigorous bound t m log n [33]. Similarly, adiabatic state preparation starting from any other state in the paramagnetic regime is efficient [24]. By contrast, the parent Hamiltonian of the zero temperature (β → ∞) Gibbs state resides at a tricritical point in the phase diagram, whose dynamical critical exponent bounds the mixing time by t m n 2 . This scaling can be understood as a consequence of the diffusive propagation of domain walls in the classical Markov chain dynamics. For the Markov chain to converge to the perfect ferromagnetic states with all spins aligned, it is necessary to remove all domain boundaries. Since the largest domains can be on the order of the system size, removing them will in general take on the order of n 2 steps in the random walk of the domain wall. The possibility of ballistic propagation in a coherent quantum system, as evidenced by the dynamical critical exponent z = 1 away from the tricritical point, suggests that adiabatic state preparation could achieve a quadratic speedup over the classical Markov chain for sampling from the Gibbs distribution of the Ising chain at zero temperature. We show in the next section that this simple argument indeed captures the relevant physics.

B. Adiabatic state preparation time
To adiabatically prepare the Gibbs state at any temperature, we start from the ground state of the Hamiltonian Eq. (12) with J 1 /h = J 2 /h = 0, corresponding to the parent Hamiltonian at infinite temperature (β = 0). The ground state is a product state of each spin aligned along the x direction, which, we assume, can be easily prepared. A measurement of this state in the computational basis yields any configuration with equal probability as required for a Gibbs distribution at infinite temperature. As argued above, adiabatic state preparation of the Gibbs sate |ψ(β) proceeds efficiently for any finite value of β. We therefore focus on the more challenging problem of sampling from the Gibbs distribution of the Ising chain at zero temperature, investigating the four adiabatic paths in Fig. 2(a). Before we proceed with the detailed analysis, we comment on an issue related to the breaking of the Z 2 symmetry by the ferromagnetic ground state. The Markov chain is in fact nonergodic at zero temperature and the gap of the generator matrix vanishes exactly. The reason for this is that both ferromagnetic configurations (all spins up or down) are stationary states of the Markov chain. It is nevertheless possible to sample from the stationary distribution by restarting the Markov chain in a random configuration. Rather than specifying a mixing time, which is ill defined, we should characterize the wait time that is necessary between restarts. Fortunately, the quantum formalism avoids this complication entirely as the system remains in the even parity subspace (see appendix A 1) at all times. The ground state of the even parity subspace remains unique, while the nonergodicity of the Markov chain is reflected in the degeneracy between the even and odd subspaces. The relevant wait time for the Markov chain can be estimated from the gap of the Hamiltonian in the even subspace. Since the distinction between the wait and the mixing time is irrelevant from the perspective of computational complexity, we will use the latter term to encompass both and we will refer to them by the same symbol t m in what follows.
To compute the adiabatic state preparation time for each path in Fig. 2(a), we numerically integrate the time-dependent Schrödinger equation, where the rate of change of the parameters is chosen following section III (see appendix A 2 for details). The resulting dependence of J 1 /h on t/t tot is plotted in Fig. 2(b) in the case of n = 100 spins. We obtain the fidelity F as a function of ε, which determines the total evolution time t tot according to Eq. (7), for different system sizes ranging from n = 10 to n = 1000. As can be seen from Fig. 3, the fidelity depends only weakly on the system size and exhibits a universal dependence on ε with 1 − F approximately proportional to ε 2 . This can be understood from the fact that the first-order diabatic correction to the wavefunction is inversely proportional to the total time t tot = l/ε [39]. From these results, we compute the adiabatic state preparation time t a as the earliest time t tot at which 1 − F is less than 10 −3 . The result is shown in Fig. 4(a).
The adiabatic state preparation time of the four paths follows distinct dependencies on system size. Path (ii), which corresponds to the one-parameter family H q (β), exhibits the scaling t a ∼ n 2 . This matches the scaling of the lower bound on the sampling time of the Markov chain at zero temperature. Paths (iii) and (iv), which cross into the ferromagnetic phase, achieve a quadratically faster scaling, t a ∼ n. This scaling is in agreement with our earlier prediction that crossing the phase transition away from the tricritical point facilitates ballistic propagation of domain walls. Finally, path (i) is slower with t a ∼ n 3 . We will see below that this slowdown can be attributed to the fact that the gap at the transition between the paramagnetic phase and the cluster-statelike phase may close exactly even for finite-sized system.
We may gain analytic insight into the above scalings by considering the adiabatic path length l. Since the fidelity is largely independent of n for fixed ε, the adiabatic path length serves as a reliable proxy for the adiabatic state preparation time t a . Indeed, Fig. 4(b) shows that l and t a exhibit roughly the same dependence on system size up to a constant prefactor. The adiabatic path length is dominated by the singular behavior close to the tricritical point. We show in appendix A 3 that in the thermodynamic limit, the adiabatic metric diverges at the tricritical point as .
Here η and α specify the distance and direction relative to the tricritical point: J 1 = 2 + η cos α, J 2 = 1 + η sin α, having set h = 1. The first case in Eq. (14) corresponds to the ferromagnetic phase (−3π/4 < α < π/4), while the second case applies to the paramagnetic and clusterstate-like phases (π/4 < α < 5π/4). In both cases, the adiabatic metric diverges as a power law, G ∼ nη −ρ , with ρ = 5/2 in the ferromagnetic phase and ρ = 5 otherwise. For finite-sized systems, one can show that exactly at the critical point, G ∼ n σ , where σ = 6 in any direction not parallel to the J 1 axis. Based on finite-sized scaling arguments, we expect that the metric follows the expression for the infinite system as we approach the tricritical point until it saturates to the final value G ∼ n σ . We are thus led to define a critical region η < η c determined by nη −ρ c ∼ n σ , in which the metric is approximately constant. These arguments imply that the path length should scale according to The above prediction agrees well with the numerical results for paths (ii)-(iv) in Fig. 4, however, it fails for path (i), for which a cubic scaling is observed. In the latter case, the scaling hypothesis breaks down since the dispersion minimum is neither at k = 0 or k = π as the tricritical point is approached from the ferromagnetic phase, thereby giving rise to an incommensurate length scale that breaks scale invariance. A similar analysis can be performed at the transition between the paramagnetic and the ferromagnetic phases, away from the tricritical point. One finds that the adiabatic path length always scales linearly with the system size, in agreement with the numerical results for paths (iii) and (iv).

A. Parent Hamiltonian
As the next example, we consider the problem of sampling from weighted independent sets. Given an undirected graph with vertices V , an independent set is any subset of V in which no pair of vertices is connected by an edge. We assign to each independent set the classical energy where the sum runs over all vertices and w i are positive weights. The occupation number n i = 1 (occupied) when the vertex i is in the independent set and n i = 0 (unoccupied) otherwise. The task at hand is to sample from the Gibbs distribution associated with H c for a given inverse temperature β. Problems of this type encompass a large class of challenging problems that appear in diverse practical settings. For instance, finding the maximum independent set, which may be viewed as sampling at zero temperature with equal weights w i = 1, has applications in computer vision [50], biochemistry [51], and social network analysis [52]. The problem is hard for general graphs even at a finite (but sufficiently low) temperature since approximately solving the maximum independent set problem on general graphs is NP hard [53,54]. The Hamiltonian in Eq. (16) is also interesting from a physics perspective as a model of a lattice gas of hard spheres [55,56].
To derive a parent Hamiltonian for the Gibbs state, we construct a Markov chain with the Metropolis-Hastings update rule [57]. A move is accepted with probability p accept = min(1, e −β∆E ), where ∆E is the change in energy. With single site updates, the probability of changing the occupation of vertex i is given by The projector P i = j∈Ni (1 − n j ) projects onto configurations where the nearest neighbors N i of vertex i are all unoccupied to ensure that the Markov chain never leaves the independent set subspace. We do not consider the possibility in which the Markov chain is initialized in a state outside this subspace. Following the same steps as in section IV A, we derive the parent Hamiltonian where V e,i (β) = e −βwi , V g,i (β) = 1, and Ω i (β) = e −βwi/2 . The Pauli operator σ x i changes the occupation number of site i. For a complete description in terms of spins, we further introduce the Pauli z operators σ z i = 1 − 2n i . We remark that constrained quantum models of this type have recently attracted great interest in the context of many-body scars [58].

B. Implementation with Rydberg atoms
Before exploring the quantum phase diagram of Eq. (18) for particular graphs, we introduce a scheme to realize this Hamiltonian in a system of trapped, neutral atoms interacting via highly excited Rydberg states. The proposal extends a previous method to encode the maximum independent set problem [59] and it is well suited for programmable atom arrays [27,28,60,61]. As a key feature of the scheme, Rydberg blockade allows for direct realization of the projectors P i , which represent d-body operators for a vertex of degree d, provided the underlying graph is a unit disk graph. This class of graphs consists of vertices embedded in two-dimensional Euclidean space, where two vertices are connected by an edge if and only if they are separated by less than a unit distance R [see Fig. 5(a)].
An atom is placed on each vertex of the graph. As illustrated in Fig. 5(b), each atom has a ground state |g i , encoding the unoccupied vertex i, and a Rydberg state |r i corresponding to the vertex being occupied. The Rydberg blockade ensures that the states of the atoms respect the independent set constraint. More concretely, we implement the first and last term in Eq. (18) by driving a transition from |g i to |r i . The value of V e,i is set by the detuning of the drive, whereas Ω i is proportional to the drive amplitude E i . The projectors P i arise naturally from Rydberg blockade: If an atom is excited to the Rydberg state, the strong van der Waals interaction U vdW shifts the Rydberg states of all neighboring atoms out of resonance, effectively turning off the drive and thereby enforcing the independent set constraint. The remaining second term in Eq. (18) can be realized using a similar approach, combining the Rydberg blockade with an AC Stark shift induced by an off-resonant drive from the ground state to an auxiliary Rydberg state |r i . The Rydberg interaction contributes additional terms to the Hamiltonian that decay as 1/r 6 with the distance r between two atoms. We have neglected these terms, noting that a strategy to mitigate their role has been proposed in a related context [62].

C. Chain graph
We now consider a chain graph of length n with periodic boundary condition and equal weights, w i = 1. The parent Hamiltonian in Eq. (18) can be written as where we dropped the subscripts i since the coefficients are translationally invariant. We further replaced P i n i by n i , which has no effect for independent sets. Equation (19) has been previously discussed as a model of strongly interacting bosons [63,64]. Since the chain graph is a unit disk graph, Eq. (19) can in principle be experimentally realized using the scheme described above. An alternative, approximate realization of the Hamiltonian has already been demonstrated by directly taking advantage of the next-to-nearest-neighbor interactions [61]. Intriguingly, an exact expression for the low-energy states can be obtained along the one-parameter family defined by β [64]. On the other hand, there is no known exact solution in the extended parameter space spanned by (Ω/V g , V e /V g ), where we assume V g > 0 throughout. We instead obtain an approximate phase diagram by numerically diagonalizing the Hamiltonian for finite chains. The complexity of the problem is reduced as we only need to consider the subspace of independent sets. Given a chain of n vertices with periodic boundary condition, the dimension of the Hilbert space is equal to F n−1 + F n+1 , where F n is the n th Fibonacci number. We can further restrict the analysis to states that are invariant under translation (zero momentum). Assuming that the initial state satisfies this condition, the state will remain in this subspace at all times since the Hamiltonian is translationally invariant. With n = 24 vertices, the dimension of this subspace is d = 4341, while for n = 30, we have d = 62075. We emphasize that the restriction to the translationally invariant subspace does not exclude spatially ordered states, as they can form a translationally invariant state by equal superpositions of the translated spatial order. It was shown in reference [63] that Eq. (19) supports distinct phases with broken translational symmetry, among them a Z 2 and a Z 3 ordered phase [ Fig. 6(b)], where the spatial order has a period of two or three vertices, respectively. To determine the approximate location of the former two phases, we introduce the order parameter for integer k. By numerically evaluating the expectation value of |M 2 | + |M 3 | in the ground state of a chain with n = 30 sites, we can clearly distinguish the Z 2 and Z 3 phases from the disordered phase [ Fig. 6(a)]. The phase boundaries are in agreement with those obtained in reference [63] with the caveat that we do not resolve the incommensurate phase that occurs in a thin slice separating the disordered phase from the Z 3 phase. The one-parameter family H q (β) is indicated by the red curve (i) in Fig. 6(a). For large β, the curve asymptotes to the phase boundary between the disordered phase and the Z 2 phase, which was determined analytically in reference [63]. The spectral gap along the oneparameter family, plotted in Fig. 6(c), allows us to bound the mixing time of the Markov chain. At high temperature, the dependence e −2β is theoretically expected [64]. Our numerical results are consistent with this prediction, although in the zero-momentum subspace it only holds for large systems. At low temperature, the gap is proportional to e −β /n 2 as can be seen from Fig. 6(c) and (d). Similar to the one-dimensional Ising model, the quadratic dependence on n is a consequence of the diffusive propagation of domain walls in the ordered phase, i.e. pairs of unoccupied sites that break up the Z 2 ordering. In contrast to the Ising model, these domain walls must overcome an energy barrier to propagate by removing the occupation of an adjacent vertex. This results in the additional factor e −β , which implies that the Markov chain is nonergodic at zero temperature even after accounting for the degeneracy that arises from spontaneous symmetry breaking. Unlike for the Ising chain, it is not sufficient to restart the Markov chain to sample from the Gibbs distribution at zero temperature.
While sampling at zero temperature is not possible, we can obtain samples that are close to the zero-temperature distribution by running the chain at a temperature that decreases with system size. The above analysis of the gap indicates that finite size effects, which are due to the correlation length reaching the system size, become relevant when e −2β ∼ e −β /n 2 . This suggests that running the Markov chain at inverse temperature β c = 2 log n will yield a high overlap with the zero-temperature Gibbs distribution. Indeed, we verified numerically that the fidelity of the ground state of the parent Hamiltonian at this temperature relative to the Gibbs state at zero temperature is approximately 90%, independent of system size. The constant overlap reflects the fact that the correlation length at β c is a fixed ratio of the system size. It is possible to increase the overlap by adding a constant to β c without changing the scaling behavior discussed in what follows. From the gap of the parent Hamiltonian, we bound the mixing time of the Markov chain at this temperature by t m n 4 .
The mixing time is again to be compared to the adiabatic state preparation time. We first observe that the Gibbs state at β = 0, which is an equal superposition of all independent sets, is not a simple product state. Nevertheless, it can be connected to a product state by a path that lies fully in the disordered phase. For example, the ground state for Ω/V g = 0, V e /V g > 3 is a product state of all sites unoccupied. We assume that this state can be readily prepared and thus, by adiabatic evolution along a path that remains in the disordered phase, any Gibbs state corresponding to a nonzero temperature (independent of n) can be prepared efficiently. The gapped spectrum of the disordered phase similarly implies that the Markov chain mixes fast for all such states.
Just as for classical sampling using the Markov chain, adiabatically preparing the zero-temperature Gibbs state  Fig. 6(a). The dashed lines are guides to the eye indicating the power laws ta ∝ n 4 and ta ∝ n. Note that the target state along path (i) is |ψ(βc) with βc = 2 log n, while path (ii) aims to reach the exact zero-temperature Gibbs state.
is challenging. The adiabatic path length l along path (i), i.e. the one-parameter family H q (β), diverges as n 3 e β/2 as shown in Fig. 7(a). It is therefore impossible to reach the zero temperature Gibbs state using the adiabatic schedule described in section III. If we follow the same strategy as for the Markov chain and instead prepare the Gibbs state at β c = 2 log n, the relevant adiabatic path length along path (i) depends on the system size as l ∼ n 4 . Intriguingly, there exist other paths, such as path (ii) in Fig. 6(a), along which the zero-temperature Gibbs state can be reached from the disordered phase with a finite adiabatic path length. Numerically, we find that the adiabatic path length along path (ii) exhibits the scaling l ∼ n [ Fig. 7(b)]. The path length is dominated by the contribution from the phase transition between the disordered phase and the Z 2 ordered phase. This transition is in the same universality class as the second order phase transition between the paramagnetic and ferromagnetic phases in the transverse-field Ising model [63]. As for the Ising model in section IV, the linear dependence of the adiabatic path length on the system size is thus a consequence of the dynamical critical exponent z = 1 or, more physically, the ballistic propagation of domain wall boundaries at the phase transition. Assuming that the adiabatic state preparation time is proportional to the adiabatic path length, as was the case for the Ising chain section IV, we expect t a ∼ n 4 along path (i) and t a ∼ n for path (ii). We verify that this assumption holds in Fig. 7(c), where we numerically integrated the Schrödinger for chains with length ranging from n = 12 to n = 24 (see appendix B 1 for details). Hence, the adiabatic state preparation time t a displayed in Fig. 7(d) exhibits the same dependence on system size as the adiabatic path length, t a ∼ n 4 for path (i) and t a ∼ n for path (ii). We reiterate that along path (i), we only aim to prepare the Gibbs state |ψ(β c ) with β c = 2 log n, which has a fidelity of approximately 90% relative to the zero-temperature Gibbs state. No such caveat applies to path (ii), for which the zero-temperature Gibbs is the target state.
To summarize, adiabatic state preparation along path (ii) achieves a quartic speedup over the classical Markov chain to sample from the Gibbs state at zero temperature. The linear scaling of the quantum algorithm with system size is explained by the ballistic propagation of domain walls, whereas the quartic scaling of the classical algorithm is the result of diffusive propagation combined with a thermal barrier at low temperature. We remark that it is possible to remove the thermal barrier in the Markov chain by introducing an update that simultaneously changes the occupation of two adjacent vertices. With such an update, one can sample from the zero-temperature Gibbs state using the Markov chain in a time t m n 2 , limited by diffusion without an energy barrier as in the Ising chain. We do not expect the pair updates to modify the scaling of the adiabatic state preparation time along path (i), which is not limited by diffusion but by the (optimal) ballistic propagation of defects. Thus, the quantum algorithm retains a quadratic speedup when pair updates are included in the Markov chain. It is intriguing that the quantum algorithm achieves the optimal scaling even though the parent Hamiltonian was derived for a suboptimal Markov chain.

D. Star graph
In the previous examples, sampling becomes hard only at zero temperature. However, sampling from the zerotemperature Gibbs distribution is equivalent to minimizing the energy of the classical Hamiltonian, for which there may exists more efficient algorithms than Markov chain Monte Carlo. It is therefore of interest to analyze a sampling problem with a Markov chain that mixes slowly at finite temperature. To this end, we consider sampling from independent sets of the star graph shown Even though the star graph is not a unit disk graph, it may be possible to implement the corresponding parent Hamiltonian using anisotropic Rydberg interactions [65]. Before discussing the Markov chain and the quantum dynamics of this model, it is helpful to consider the classical equilibrium phase diagram. The partition function corresponding to the classical Hamiltonian H c in Eq. (16) is given by Z = 1 + 2e β b +e bβ 1 + e β b . The two terms arise from the different configurations of the central vertex. From the Helmholtz free energy F = − log Z/β and the total energy U = −∂ log Z/∂β, the entropy can be computed as S = β(U − F ), which is plotted in Fig. 8(c). In the thermodynamic limit, the system undergoes a discontinuous phase transition at β * = log ϕ, where ϕ = ( √ 5 + 1)/2 is the golden ratio. The origin of the phase transition can be understood by noting that the probability that the central site is occupied is given by This expression turns into the step function p 1 = Θ(β − β * ) in the thermodynamic limit b → ∞. At high tem-perature, it is entropically favorable for the central vertex to be unoccupied as this allows the branches to be in 3 b distinct configurations. Below the phase transition temperature, the reduction in energy from occupying the central vertex outweighs the entropic cost of reducing the number of configurations to 2 b due to the independent set constraint. The Markov chain with single spin flips on this graph is subject to a severe kinetic obstruction since changing the central vertex from unoccupied to occupied requires all neighboring vertices to be unoccupied. Assuming that each individual branch is in thermal equilibrium, the probability of accepting such a move is given by p 0→1 = [(1 + e β )/(1 + 2e β )] b . The reverse process is energetically suppressed with an acceptance probability p 1→0 = e −bβ . The central vertex can thus become trapped in the thermodynamically unfavorable configuration, resulting in a mixing time that grows exponentially with the number of branches b at any finite temperature. Indeed, we will see below that the spectral gap of the parent Hamiltonian (and thus the Markov chain) is approximately equal to p 1→0 for temperatures above the phase transition and to p 0→1 below. We remark that despite the exponentially large mixing time, it is possible to sample efficiently at high temperature if the initial configuration for the Markov chain is drawn from the uniform distribution. Since the probability of the central vertex being initially occupied is exponentially small, the desired Gibbs distribution can be approximated with exponentially small error by restarting the Markov chain multiple times. By the same argument, the Markov chain almost certainly starts in the wrong configuration in the low temperature phase and the exponentially large mixing time t m 1/p 0→1 is a limiting factor.
To determine the quantum phase diagram associated with the parent Hamiltonian Eq. (18) for the star graph, we observe that the parent Hamiltonian is invariant under permutations of the branches. Restricting the analysis to the completely symmetric subspace, we show in appendix C that H q (β) possesses a two-dimensional lowenergy subspace. The subspace is to an excellent approximation spanned by the states |ψ 0 (β) and |ψ 1 (β) , which represent the Gibbs states with the central vertex respectively fixed to be unoccupied or occupied. This is true even when the parameters Ω i and V e,i associated with the central vertex, denoted by Ω cen and V e,cen , are varied arbitrarily, assuming all other terms in the Hamiltonian follow the one-parameter family H q (β). The effective Hamiltonian in the low-energy subspace can be obtained by a Schrieffer-Wolff transformation [66], which to second order yields where The term f is a second order correction, which vanishes as 1/b in the thermodynamic limit.
Since |ψ 0 (β) and |ψ 1 (β) represent macroscopically distinct states, the system undergoes a first-order quantum phase transition when the bias (1 − f )ε 0 − (V e,cen − f Ω 2 cen ) vanishes. Along the one-parameter family H q (β), we have V e,cen = Ω 2 cen such that the phase transition occurs when ε 0 = V e,cen . Solving for β gives the critical value β * = log ϕ, in agreement with the location of the classical phase transition. More generally, Eq. (22) predicts a gap given by ( As shown in Fig. 8(d), this expression agrees well with the numerically exact result of the gap in the permutation symmetric subspace for large values b. The two-state model thus tells us that the gap vanishes as e −bβ above the phase transition (β < β * ) and as (1 + e β ) b /(1 + 2e β ) b below (β > β * ) up to a small correction due to the (1−f ) prefactor. This result perfectly matches the estimates of the mixing time of the Markov chain based on the transition probabilities p 0→1 and p 1→0 . It is worth pointing out that the location of the quantum phase transition is well defined even though the gap vanishes for any β > 0 in the thermodynamic limit. We can define the critical region of the phase transition as the range of parameters for which the tunneling rate (1 − f )J exceeds the bias The critical region is of size ∆β ∼ 1/b, such that the phase transition point β * is indeed meaningful.
Having described the Markov chain and the related spectral properties of the parent Hamiltonian H q (β), we proceed to the adiabatic state preparation of the Gibbs state. We note that the Gibbs state at β = 0 can be efficiently prepared from a product state. For instance, one can start with V e,i = 1 and V g,i = Ω i = 0, where the ground state is the state of all vertices unoccupied. Next, all Ω i are simultaneously ramped up to 1 at a constant rate before doing the same for V g,i . The Hamiltonian remains gapped along this path such that the adiabatic state preparation proceeds efficiently. Hence, we may use the Gibbs state at β = 0 as the initial point for all adiabatic paths. Despite the exponentially small gap, it is also possible to efficiently prepare any state above the phase transition. The probability of occupying the central vertex is exponentially small above the critical temperature, outside the critical region. Therefore, an exponentially good approximation to the Gibbs state can be prepared along the one-parameter family H q (β) since adjusting the amplitudes of the different configurations of the branches proceeds without kinetic obstruction.
As for the Markov chain, this line of reasoning fails when the target temperature is below the phase transition. Along the one-parameter family H q (β), the gap at the phase transition point is equal to the tunneling rate J = ϕ −b , ignoring the subleading factor (1 − f ). Hence, we expect the adiabatic state preparation time to be proportional to ϕ b , which we numerically confirmed for the final state |ψ(2β * ) in Fig. 9. The numerical integration of the Schrödinger equation follows the same steps as for the independent set problem on the chain graph (see appendix B 1) with time steps chosen to satisfy ∆t||d|0 /dt|| = 10 −4 . We note that the lower bound on the mixing time of the Markov chain at the phase transition has the same scaling t m ϕ b . For sampling below the phase transition, the mixing time increases such that the quantum algorithm appears to achieve a speedup. However, the slowdown of the Markov chain can be circumvented by performing simulated annealing [2] across the phase transition. While a detailed analysis of simulated annealing is beyond the scope of the present work, we expect that its convergence time matches the time to adiabatically prepare a state along the one-parameter family H q (β) for any final value of β since the occupation of the central vertex is essentially frozen outside the critical region of the phase transition. The adiabatic state preparation time along the oneparameter family H q (β) is limited by the tunneling rate J at the phase transition. Since J is proportional to Ω cen , it is natural to consider trajectories along which Ω cen is of order one at the phase transition rather than e −bβ * /2 . Because [(1+e β * )/(1+2e β * )] b = e −bβ * , this simple argument together with Eq. (23) predicts a quadratic speedup. As a first guess, one may consider a path where Ω cen is held constant at 1, while all other parameters are varied according to H q (β) from β = 0 to the desired final value of β, taken to be 2β * in what follows. The time required to cross the phase transition is again t a ∼ 1/J, where J is evaluated at the phase transition. The term f Ω 2 cen in the effective Hamiltonian Eq. (22) shifts the phase transition close to β = 0 such that t a ∼ (3/2) b/2 . However, this does not yet result in a speedup for preparing the Gibbs state as one still has to decrease Ω cen to its final value e −bβ * . It turns out that this step negates the speedup if performed adiabatically. The reason for this is that the large value of Ω cen admixes states where the central vertex is unoccupied, skewing the occupation probability away from its thermal expectation value. Even though this admixture is small, on the order of f , the time to adiabatically remove it exceeds the time to initially cross the phase transition. To avoid this issue, one can in principle suddenly switch Ω cen to its final value at the cost that the final fidelity is limited by the admixture. Perturbation theory and numerical results indicate that this results in an infidelity 1 − F that decays as 1/b 2 for large values of b.
A slight modification of the path achieves a final state infidelity that is not only polynomially but exponentially small in b. First, V e,cen is lowered from its initial value 1 to −1, which can be done in time t a ∼ (3/2) b/2 as before. Next, all other parameters are varied along H q (β), for which only a time polynomial in b is required. Numerical results for these two steps are shown in Fig. 10. Finally, V e,cen is ramped to its final value e −2bβ * . Similar to the previous scheme, perfectly adiabatic evolution of this last step would require a very long time because the probability that the central vertex is unoccupied is initially too small due to the large, negative value of V e,cen . However, since its final value p 0 ≈ e −2bβ * [(1 + 2e 2β * )/(1 + e 2β * )] b is exponentially small in b, the fidelity can be exponentially close to unity when suddenly switching V e,cen to its final value. We numerically confirmed this prediction for the sudden transition in Fig. 11.
To summarize, we identified a suitable path along which the adiabatic evolution, supplemented by a sudden change of parameters at the end, achieves a preparation time of t a ∼ (3/2) b/2 . The speedup over the Markov chain with mixing time t m ϕ b at the phase transition appears to be more than quadratic. However, it is likely that a classical computation time on the order of (3/2) b can be achieved using simulated annealing. For instance, one could consider an annealing schedule in which the weight on the central vertex is first increased. This shifts  Fig. 10. The value of Ve,cen is switched suddenly from −1 to e −2bβ * . The dots show the resulting infidelity, assuming perfect fidelity along the preceding adiabatic path. The infidelity agrees well with the probability that the central vertex is unoccupied (dashed line, whose functional form is shown in the plot).
the phase transition towards β = 0, allowing one to sample at the phase transition in a time that scales as (3/2) b . In the annealing schedule, the temperature can then be lowered to the desired value before ramping the weight of the central vertex back to its initial value. This annealing schedule is in many ways similar to the adiabatic path discussed above. It is nevertheless quadratically slower because unlike in the quantum case it is not possible to vary V e,cen and Ω cen independently.
The above results for the star graph exhibit important differences from the chain graph and from sampling from the Ising model. While the speedup is also quadratic, it originates from coherent tunneling as opposed to the ballistic propagation of domain walls. This is related to the fact that the parent Hamiltonian for the star graph exhibits a first-order quantum phase transition as opposed to the second-order transition in the preceding models. In the next section, we discuss another model in which the dynamics are governed by a first-order transition.

VI. UNSTRUCTURED SEARCH PROBLEM
The unstructured search problem was pivotal in the development of quantum algorithms. Grover's algorithm gave an early example of a provable quantum speedup and it remains an essential subroutine in many proposed quantum algorithms [67]. Moreover, the unstructured search problem played a crucial role in the conception of adiabatic quantum computing [68]. Below, we show that when applied to the unstructured search problem, our formalism recovers the adiabatic quantum search algorithm along with its quadratic speedup over any classical algorithm. While the nonlocality of the resulting parent Hamiltonian render it challenging to implement in practice, the result underlines the power of our approach in enabling quantum speedup in a setting where the speedup is provably optimal [69].
We consider the problem of identifying a single marked configuration m in a space of a total of N elements. To connect this search problem to a sampling problem, we assign the energy −1 to the marked configuration, while all other states have energy 0. This is summarized by the classical Hamiltonian Solving the search problem may now be formulated as sampling from the Gibbs distribution associated with H c at zero temperature. Given the lack of structure of the problem, a natural choice for the Markov chain is to propose any configuration with equal probability 1/N [70]. If the update is accepted according to the Metropolis-Hastings rule, the resulting parent Hamiltonian takes the form where Since |ψ 0 is contained in the subspace spanned by {|m , |m ⊥ }, the Hamiltonian acts trivially on the orthogonal subspace. In fact, all nontrivial dynamics arise from the second line of Eq. (25). Considering the extended parameter space of arbitrary (V 0 , V m ), we can show that the Hamiltonian supports two distinct quantum phases separated by a first-order phase transition. By rewriting Eq. (25) in terms of |m and |m ⊥ only, one can show that these two states have a relative energy difference V m − (1 − 2/N )V 0 while being connected by an off-diagonal matrix element (tunneling rate) of magnitude √ N − 1V 0 /N . In the thermodynamic limit N → ∞, the system thus undergoes a first-order quantum phase transition when V m = V 0 . For V m > V 0 , the ground state has large overlap with |m , whereas for V m < V 0 , it is close to |m ⊥ (and |ψ 0 ).
The one-parameter family H q (β), indicated by the red curve (i) in Fig. 12(a), traces out a segment of a parabola passing through (V 0 , V m ) = (1, 0) when β = 0 and (V 0 , V m ) = (0, 1/N ) as β → ∞. The gap at zero temperature is equal to 1/N [ Fig. 12(b)], which allows us to bound the mixing time by t m N (up to logarithmic corrections). This bound is expected as any classical algorithm must check on average half the configurations to solve the unstructured search problem. As in all previous examples, adiabatic state preparation along the one-parameter family leads to the same time complexity. To see this, we assume that the ground state at β = 0, given by |ψ 0 , can be readily prepared. Adiabatic state preparation experiences a bottleneck close to the quantum phase transition, where V 0 and V m are on Since the adiabatic state preparation time is limited by the inverse of the tunneling rate in the critical region, we obtain an adiabatic state preparation time t a ∼ √ N /V 0 ∼ N . Similar to what we found for the star graph, the adiabatic state preparation can be sped up by crossing the phase transition at a point where the tunneling rate is large. In particular, this requires that V 0 and V m be of order one. One such path is path (ii) shown in Fig. 12(a). A straight line segment connects (V 0 , V m ) = (1, 0) to (0, 1) before continuing to the final point (V 0 , V m ) = (0, 1/N ). Since the Hamiltonian is purely diagonal along V 0 = 0 in the computational basis, there are no diabatic transitions along the latter segment and the parameters can be changed suddenly. The gap along the former segment is shown in Fig. 12(b). The corresponding Hamiltonian is in fact identical to the Hamiltonian of the adiabatic quantum search algorithm, which was derived by interpolating between the projectors |ψ 0 ψ 0 | and |m m| [68]. It was shown that by carefully choosing the rate of change using a scheme essentially equivalent to that outlined in section III, it is indeed possible to prepare the final ground state with high fidelity in a time t a ∼ √ N limited by the tunneling rate at the phase transition.

VII. SUMMARY AND OUTLOOK
In summary, we have described a method to construct quantum algorithms to sample from Gibbs distributions. The approach can be readily generalized to any probability distribution that can be described as the stationary distribution of a Markov chain satisfying detailed balance. Our results differ from previous work by considering adiabatic state preparation in a parameter space that has been extended beyond the one-parameter family of Hamiltonians H q (β). By means of four examples, we showed that it is possible to achieve a quantum speedup by suitably navigating the quantum phases in the extended parameter space. The speedup has a different origin depending on the nature of the phase transition. In the case of second-order phase transitions, the speedup was due to the ballistic propagation of domain walls as opposed to diffusive motion in the classical Markov chain. For first-order phase transitions, we could trace the speedup to coherent tunneling between macroscopically distinct states.
The quantum Hamiltonians encountered in our construction are guaranteed to be local provided that the Gibbs distribution originates from a local classical Hamiltonian and that the Markov chain updates are local. This was the case in all of the examples except for the unstructured search problem, which we included to highlight the power of the approach in a more abstract setting. It is therefore possible to efficiently implement time evolution under these Hamiltonians on a universal quantum computer using Hamiltonian simulation [71]. Moreover, for sampling from independent sets of unit disk graphs, there exists a natural implementation of the parent Hamiltonian using Rydberg states of neutral atoms. The proposed scheme is compatible with existing architectures [27,28,60,61], opening the door to exploration of sampling problems on near-term quantum devices.
Further work is required to extend the applicability of our approach to a wider range problems. As a first step, one may consider a generalization of the above spin models to higher dimensions. For instance, Glauber dynamics in the two-dimensional Ising model differs substantially from its one-dimensional counterpart owing to presence of a finite-temperature phase transition in the classical model. At temperatures above the phase transition, the Markov chain still mixes rapidly, whereas at low temperature, the mixing time diverges exponentially with the linear dimension of the system [72]. Moreover, there exist many configurations with large domains, which relax slowly to equilibrium in the Markov chain. Below the phase transition, the corresponding parent Hamiltonian H q (β) therefore describes a a system with a large number of states energetically close to the ground state, hinting at the presence of an unusual quantum phase. Additional research is needed to fully characterize this quantum phase, e.g. using quantum Monte Carlo techniques, and to identify possible mechanisms for quantum speedup when adiabatically approaching it.
Practically relevant, hard sampling problems, such as sampling from the Gibbs distribution of a classical spin glass or other disordered models in two or more dimensions, lack much of the structure such as translational symmetry of the problems discussed above. Here, a key challenge is to identify suitable adiabatic paths as it may not be possible to determine the complete phase diagram. This issue could potentially be addressed by employing hybrid algorithms which combine quantum evolution with classical optimization to identify a good adiabatic path [73]. More concretely, the energy of the parent Hamiltonian H q (β) could be minimized using varia-tional quantum algorithms similar to existing proposals but without the need for complex measurements of the entanglement entropy [74,75]. Since variational algorithms do not require a physical realization of the parent Hamiltonian, this approach could be particularly fruitful for complex cost functions composed of several parent Hamiltonians involving nonlocal cluster updates such as those of the Swendsen-Wang algorithm [76]. The Hamiltonian in Eq. (12) can be mapped onto a free-fermion model using a Jordan-Wigner transformation. We define the fermion annihilation and creation operators a i , a † i and relate them to the Pauli matrices according to Equation (12) becomes, up to a constant, where N = n i=1 a † i a i is the total number of fermions. While the fermion number itself is not conserved, the parity e iπN is, allowing us to consider the even and odd subspaces independently.
We define the momentum space operators which satisfy fermionic commutation relations for suitably chosen k. We let for l = 0, 1, . . . , n − 1 (mod n). With this definition, the inverse Fourier transformed operators have the formal property a i+n = −e iπN a i , which accounts for the boundary terms in Eq. (A4). The Hamiltonian simplifies to While the above Hamiltonian can be diagonalized by a standard Bogoliubov transformation, it will prove more convenient for our purposes to map it onto noninteracting spins. For 0 < k < π, we define It is straightforward to check that these operators satisfy the same commutation relations as Pauli matrices. In addition, operators corresponding to different values of k commute such that we can view them as independent spin-1/2 systems, one for each value of k. We restrict the range of momenta to 0 < k < π due to the redundancy τ α −k = −τ α k . The cases k = 0 and k = π require special treatment as both τ x k and τ y k vanish. For concreteness, we assume that the number of spins n is even. The special cases k = 0 and k = π are then both part of the odd parity subspace (e iπN = −1). The Hamiltonian of the even parity subspace can be written as where The angles θ k are uniquely defined by The ground state is given by where |vac is the vacuum with respect to the a k operators. The ground state energy is In the odd parity subspace, we have The construction of the ground state is analogous to the even case with the additional requirement that either the a 0 fermion or the a π fermion, whichever has the lower energy, be occupied. One can show that the resulting energy is gapped above E even 0 when h + J 1 + J 2 and h − J 1 + J 2 have the same sign. In the case of opposite signs, the even and odd sector ground states are degenerate in the thermodynamic limit, corresponding to the symmetry breaking ground states of the ferromagnetic phase.
In the main text, we consider adiabatic evolution starting from the ground state at J 1 = J 2 = 0. Following the above discussion, this state is part of the even subspace. Since the time evolution preserves parity, we may restrict our discussion to the even subspace, dropping all associated labels in what follows.
Excited states can be constructed by flipping any of the τ spins. Since any spin rotation commutes with the parity operator, singly excited states are given by with an energy 4E k above the ground state.

Time-dependent Schrödinger equation
To compute the fidelity, we numerically integrate the Schrödinger equation for each spin τ k . We work in the instantaneous eigenbasis |χ ± k (t) , which are eigenstates of H k = 2E k (cos θ k τ z k + sin θ k τ y k ) with energies ±2E k [see Eq. (A11)]. It is convenient to parametrize each adiabatic path by a dimensionless time s running from 0 to 1. Writing the state at time s as the coefficient c k and d k are determined by the Schrödinger equation with the initial condition c k (0) = 1, d k (0) = 0. The final fidelity is obtained by solving this equation for each spin and multiplying the individual fidelities, We note that all terms in Eq. (A20) can be evaluated without having to solve for the physical evolution time t(s). The terms E k (s) and dθ k /ds are readily computed from equations (A12)-(A14), while dt/ds follows from equation (5) of the main text: Here, λ 1 = J 1 , λ 2 = J 2 , setting h = 1 throughout. To vary the total evolution time t tot , we simply adjust the value of ε. We obtained good convergence by evolving under constant s = s n for an interval ∆s n = 2×10 −3 / dθ k ds s=sn before incrementing s n+1 = s n +∆s n . The number of steps is independent of the total time, yet the final fidelity is well estimated since the probability of leaving the ground state is small in each step.
In the thermodynamic limit, the momentum sum turns into an integral, which can be evaluated analytically. By expanding around the tricritical point, we obtain Eq. (14). Table I gives an explicit parametrization of the paths (i)-(iv) in Fig. 2(a). The parameter s ranges from 0 to 1. For path (ii), s is related to β by s = tanh β.

Adiabatic paths
Appendix B: Sampling from weighted independent sets

Numerical details
To integrate the Schrödinger equation, we exactly diagonalize the Hamiltonian at discrete time steps ∆t. 2s s 2 (iii) 3(1 − s) 2 s + 7.5(1 − s)s 2 1.5(1 − s)s 2 + s 3 +2s 3 (iv) 6(1 − s) 2 s + 9(1 − s)s 2 −3(1 − s) 2 s + 4.5(1 − s)s 2 +2s 3 +s 3 The time steps are chosen such that ∆t||d|0 /dt|| = 10 −3 , where |0 denotes the instantaneous ground state. The expression is most conveniently evaluated using the identity ||d|0 /dt|| 2 = n>0 | n|dH/dt|0 | 2 /(E n − E 0 ) 2 , where the sum runs over all excited states |n with energy E n . The time is related to the parameters of the Hamiltonian by Eq. (5). Table I gives an explicit parametrization of the paths (i) and (ii) in Fig. 6(a). The parameter s ranges from 0 to 1. Along path (i), s and β are related by s = e −β/2 . The star graph has three types of vertices: the vertex at the center and the inner and outer vertices on each branch. If we maintain the permutation symmetry between the branches, the parent Hamiltonian takes the general form H q = V e,cen n cen + V g,cen P cen (1 − n cen ) − Ω cen P cen σ x

Adiabatic paths
where each row relates to a separate type of vertex and the sums run over all branches. With the weights specified in the main text, we have V e,cen = e −bβ , V e,in = V e,out = e −β , V g,cen = V g,in = V g,out = 1, Ω cen = e −bβ/2 , Ω in = Ω out = e −β/2 along the one-parameter family H q (β). We restrict our analysis to the subspace that is completely symmetric under permutations of the branches. We introduce the total occupation numbers n in = b i=1 n in,i and n out = b i=1 n out,i as well as the number of unoccupied branches n 0 . The symmetric subspace is spanned by the states |n cen , n in , n out , n 0 , where n cen ∈ {0, 1} while the other occupation numbers are nonnegative integers satisfying n in + n out + n 0 = b. If n cen = 1, the independent set constraint further requires n in = 0. Each of the states in Eq. (C2) is an equal superposition of b!/(n in ! n out ! n 0 !) independent configurations. The dimension of the completely symmetric subspace is (b + 1)(b + 4)/2. The permutation symmetry leads to a bosonic algebra. We define the bosonic annihilation operators b in , b out , and b 0 , respectively associated with the occupation numbers n in , n out , and n 0 , which may be viewed as a generalization of Schwinger bosons. We split up the Hamiltonian into blocks where the central spin is either 0 or 1 as well as an off-diagonal term coupling them. Explicitly, where P (n in = 0) projects onto states with no occupied inner vertices. We diagonalize the Hamiltonian by treating the projectors perturbatively. We focus on the situation where all parameters follow the one-parameter family H q (β) except for Ω cen and V e,cen , which may be adjusted freely. By diagonalizing the matrices in Eq. (C4) and Eq. (C5), we identify the lowest energy modes of the quadratic parts of H (0) q and H (1) q and associate with them the bosonic annihilation operators c 0 and c 1 , respectively. Both modes have zero energy while the other modes are gapped at any finite value of β. We may thus expect the ground state to be well approximated in the subspace spanned by where |vac denotes the bosonic vacuum. One can show that these states correspond to the Gibbs state of the star with the central spin held fixed. Next, we perform a Schrieffer-Wolff transformation [66] to project onto the subspace spanned by |ψ 0 and |ψ 1 . We arrive at an effective Hamiltonian where the terms where we used the relation P (n in = 0)|ψ 0 = √ ε 0 σ x cen |ψ 1 , which holds along the paths of interest. The sums run over all excited states |n with energy E n of the unperturbed part of H (0) q . We neglected a term V e,cen in the energy denominator, which is justified as long as V e,cen is small compared to E n . The discussion remains valid even if this is not the case because the shifts from the Schrieffer-Wolff transformation can then be ignored as far as the ground state is concerned.
The complete effective Hamiltonian may be written as where f = n | n|σ x cen |ψ 1 | 2 /E n . We find numerically that f decays as an inverse power law in b such that our approximations are well justified in the thermodynamic limit. Along the one-parameter family H q (β), we have V e,cen = Ω 2 cen such that H eff depends on f only through an overall factor (1 − f ), which tends to 1 in the limit of large b.