Introduction

Adiabatic quantum computing (AQC) is an alternative approach to traditional gate-based quantum computing where quantum adiabatic evolution is performed in order to achieve a computation1,2,3,4. In the scheme, the aim is to find the ground state of a Hamiltonian HZ which encodes the problem to be solved and can be considered an instance of quantum annealing5,6,7,8. In addition, an initial Hamiltonian HX, which does not commute with the problem Hamiltonian, is prepared such that the ground state is known. For example, a common choice of these Hamiltonians are

$${H}_{Z}^{\,\text{(qubit)}\,}=\mathop{\sum }\limits_{i=1}^{M}\mathop{\sum }\limits_{j=1}^{M}{J}_{ij}{\sigma }_{i}^{z}{\sigma }_{j}^{z}+\mathop{\sum }\limits_{i=1}^{M}{K}_{i}{\sigma }_{i}^{z},$$
(1)
$${H}_{X}^{\,\text{(qubit)}\,}=-\mathop{\sum }\limits_{i=1}^{M}{\sigma }_{i}^{x},$$
(2)

where \({\sigma }_{i}^{x,z}\) are Pauli matrices on site i, and Jij and Ki are coefficients which determine the problem to be solved, and there are M qubits. The form of Eq. (1) directly encodes a wide variety of optimization problems, for example, MAX-2-SAT and MAXCUT which are NP-complete problems. It then follows that any other problem in NP can be mapped to it in polynomial time9,10,11. AQC then proceeds by preparing the initial state of the quantum computer in the ground state of HX, then applying the time-varying Hamiltonian

$$H=(1-\lambda (t)){H}_{X}+\lambda (t){H}_{Z},$$
(3)

where λ(t) is a time-varying parameter that is swept from 0 to 1. Intense investigation into the performance of AQC has been performed since its original introduction, demonstrating its performance for various problems12,13,14,15 and characterizing the effects of decoherence16,17,18,19,20,21,22.

In the AQC framework, the speed of the computation is given by how fast the adiabatic sweep is performed. To maintain adiabaticity, one must perform the sweep sufficiently slowly, such that the system remains in the ground state throughout the evolution. The sweep time required to maintain adiabaticity is known to be proportional to a negative power of the minimum energy gap of the Hamiltonian (Eq. 3), where the power depends upon the annealing schedule λ(t) and gap structure2,13,15,23,24,25,26. One of the attractive features of AQC is that time-sequenced gates do not need to be applied, but it is nevertheless known to be equivalent to the gate-based quantum computation13,27,28,29,30,31,32. Numerous theoretical analyses and experimental demonstrations have been performed both at small33,34,35,36,37 and larger scale38,39,40,41,42,43,44,45,46. One of the outstanding problems for AQC is to fully understand the performance and effectively combat decoherence in AQC such that it can be applied to real-world combinatorial problems4,47,48.

In this paper, we investigate a variant of the AQC Hamiltonians (Eq. 1) and (Eq. 2) where ensembles of qubits are used to encode the optimization problem, instead of genuine qubits. Specifically, we study the Hamiltonians

$${H}_{Z}=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{M}\mathop{\sum }\limits_{j=1}^{M}{J}_{ij}{S}_{i}^{z}{S}_{j}^{z}+\mathop{\sum }\limits_{i=1}^{M}{K}_{i}{S}_{i}^{z},$$
(4)
$${H}_{X}=-\mathop{\sum }\limits_{i=1}^{M}{S}_{i}^{x},$$
(5)

where \({S}_{i}^{x,z}=\mathop{\sum }\nolimits_{n = 1}^{N}{\sigma }_{i,n}^{x,z}\) are total spin operators for an ensemble consisting of N qubits. In this case, M is the number of ensembles. Here, σi,n denotes the Pauli operator for the nth qubit within the ith ensemble. The matrices Jij and Ki are the same as that in Eq. (1) and we take Jij = Jji and Jii = 0. AQC then proceeds in the same way as described in Eq. (3). Each of the ensembles is initially prepared in a fully polarized state of \(\langle {S}_{i}^{x}\rangle =N\) and adiabatically evolved to the ground state of HZ. The aim will be to investigate whether the ensemble version of the Hamiltonian can be used in place of the qubit Hamiltonian, such that the ground state configuration of Eq. (1) is found using Eqs. (4) and (5). We characterize the nature of the ground and excited states of the ensemble Hamiltonian and assess the performance of AQC in comparison to the original qubit problem.

The Hamiltonians (Eq. 4) and (Eq. 5) can be considered an error-suppressing encoding of the original AQC Hamiltonians (Eq. 1) and (Eq. 2), respectively. The use of an ensemble duplicates the quantum information since the N qubits within an ensemble are in the same state at the start and at the end of the adiabatic evolution. Such error-suppression strategies have been already investigated in the context of AQC. Jordan, Farhi, and Shor49 introduced an encoding capable of detecting the presence of single-qubit errors, which are suppressed by an additional energy penalty term in the total Hamiltonian. Pudenz, Lidar, and co-workers50 introduced a scheme known as quantum annealing correction (QAC), where a repetition code is used to encode a logical qubit and a majority vote is used to decode the logical operations. This, combined with the addition of energy penalty terms to enforce alignment within each ensemble has shown that the scheme is effective at achieving error-suppression51,52,53,54,55. In particular, the two-qubit interaction term in Eq. (1) is encoded by repeating the interactions N times in a pair-wise fashion between the logical qubit ensembles, motivated by the chimera qubit connectivity of the D-wave quantum computers. An alternative scheme called nested QAC has also been studied by Vinci, Lidar, and co-workers, where the qubit terms are mapped with an all-to-all interaction, similar to the one we consider in Eq. (4)56. A similar all-to-all encoding was considered by Venturelli, Smelyanskiy, and co-workers57. A minor embedding is then performed on the encoded Hamiltonian to match it to the chimera graph topology58. Matsuura and co-workers showed through mean-field analysis that the order of the phase transition is modified through the mapping in the ferromagnetic and antiferromagnetic Ising models59. Young, Sarovar, and Blume-Kohout47,48 showed that such energy penalty approaches could be part of a unified theory with dynamical decoupling.

Results

Properties of the problem Hamiltonian H Z

We first examine the properties of the ensemble version of the (classical) problem Hamiltonian (Eq. 4). The typical energy landscape of the Hamiltonian HZ is shown in Fig. 1a. The axes are plotted with rescaled spin variables \({x}_{i}=\langle {S}_{i}^{z}\rangle /N\). For the corners of the hypercube xi = ± 1, the energy eigenvalues of the ensemble Hamiltonian HZ reduces to that of the qubit version (Eq. 1) up to an overall scaling factor of N. This is a general result which is true by virtue of the structure of the ensemble and qubit problem Hamiltonians HZ (see Supplementary Note 1). The primary difference between Eq. (4) and Eq. (1) is then that the ensemble version can take a discrete set of intermediate values of xi between the ±1 values.

Fig. 1: Properties of the problem Hamiltonian HZ.
figure 1

a Summary of the two regimes for the ensemble version of the problem Hamiltonian HZ. The cube shows the energy landscape of a typical instance of Eq. (4) with rescaled variables \({x}_{i}=\langle {S}_{i}^{z}\rangle /N\). Dots indicate allowed states for N = 1 (left) and N = 6 (right). Only states on the visible faces are shown for clarity. b Energy variation for the problem instance as Fig. 1a along a linear trajectory from the ground state x1 = x2 = x3 = 1 to other hypercube corners specified by \(({x}_{1}^{(f)},{x}_{2}^{(f)},{x}_{3}^{(f)})\). The trajectory is defined by \({x}_{i}=1-(1-{x}_{i}^{(f)})\varepsilon\), where ε = 0 is the ground state and ε = 1 is another hypercube corner. Parameters used are M = 3, J12 = − 2, J13 = − 1, J23 = − 1, K1 = 1, K2 = 0, K3 = − 2. c Proportion of 800 randomly generated problem instances with Δ < δ (i.e., N < Nc) as a function of the ensemble size N. Elements of Jij and Ki are taken randomly in the interval \(\left[-1,1\right]\) with a uniform distribution.

Full size image

In Fig. 1b, we show the variation of the energy starting from the ground state to the remaining hypercube corners. The variation always follows a quadratic form with an initially positive slope. This can be shown to be generally true following from the quadratic form of Eq. (4) (see Supplementary Note 1). This fact can be used to show that any point along a trajectory connecting the ground state to another hypercube corner has energies greater than the ground state (see Supplementary Note 1). From the above structure of the energy landscape, one can deduce that the ground states of the Hamiltonians (Eq. 1) and (Eq. 4) have logically equivalent spin configurations. We define the logically equivalent states of the qubit and ensemble systems according to

$$\,\text{sgn}\,[{\langle {\sigma }_{i}^{z}\rangle }^{\text{(qubit)}}]=\,\text{sgn}\,[{\langle {S}_{i}^{z}\rangle }^{\text{(ens)}}],$$
(6)

for all i. This is also known as a majority vote encoding of the ensemble to give the logical state and has been considered in other error mitigation schemes for AQC47,50,56. Thus, finding the ground state spin configuration of Eq. (4) gives the same information as Eq. (1).

In AQC, one of the parameters which play a central role is the gap energy, i.e., the energy between the ground and first excited state. The simple structure of the Hamiltonian (Eq. 4) allows us to deduce that for a given N, there are two different regimes for the gap energy. For the particular example shown in Fig. 1a), we see that there are two hypercube corner states (x1,x2,x3) = (1,1,1) and (−1,−1,−1) with relatively similar energies of ϵ0 and ϵ1, respectively. For the qubit case (Eq. 1) the difference between the two lowest energy hypercube corners is the gap energy. In the ensemble case, the energy difference between these two hypercube corners is

$${{\Delta }}=N({\epsilon }_{1}-{\epsilon }_{0}),$$
(7)

since for extremal values ∣xi∣ = 1 the energies are related by a factor of N.

If N is sufficiently small, Δ remains the gap energy for the ensemble case. However, for large enough N this becomes less and less likely, and the first excited state is a single qubit spin-flip of the ground state. Specifically, we have the state such that on the kth ensemble \({S}_{k}^{z}=\pm (N-2)\), and the remaining ensembles \({S}_{i}^{z}=\pm N,\forall i\ne k\). The energy gap for this single qubit flip state is

$$\delta =\mathop{\min }\limits_{k}\left[-2{\sigma }_{k}\left({K}_{k}+2\mathop{\sum}\limits_{j}{J}_{kj}{\sigma }_{j}\right)\right],$$
(8)

where \({\sigma }_{i}={\langle {\sigma }_{i}^{z}\rangle }_{\,\text{gnd}}^{\text{(qubit)}}=\text{sgn}\,({\langle {S}_{i}^{z}\rangle }_{\,\text{gnd}}^{\text{(ens)}\,})\), and expectation values are taken with respect to the ground state. The minimum function picks the smallest value from the range k ∈ [1, M].

Whether Δ or δ is the gap energy depends upon N and the particular parameter choice of Jij and Ki made. As N is increased, at some point there will always be a crossover such that the gap is δ, since (Eq. 7) is proportional to N while (Eq. 8) has no dependence on N. Let us call Nc the smallest value of N such that Δ > δ. For a given problem instance, we then define two regions of N, according to whether it is larger or smaller than Nc. The two regimes and their implications summarized in Fig. 1a. In Fig. 1c we show the proportion of problems satisfying Δ < δ (i.e., N < Nc) for randomly generated Jij and Ki. We see that the proportion decreases as ∝1/N, which is consistent with the linear scaling of Δ. We note that in the case of atomic ensembles N can be quite large (e.g., 103 to 1011)60,61, which suggests that for realistic ensemble sizes most of the problem instances will be in the regime N≥Nc.

The spectrum of the AQC Hamiltonian

So far we have only examined the classical limit of λ = 1. The overall speed of the AQC will be dependent upon the minimum gap energy with the off-diagonal term (Eq. 5) present. To illustrate the effect of intermediate λ, we compare the eigenvalue spectrum of the Hamiltonian (Eq. 3) for the standard qubit case and the ensemble case with N = 5 for the same Jij and Ki parameters in Fig. 2a, b. Due to symmetry under particle interchange on each ensemble, the Hilbert space can be reduced to the symmetric subspace, reducing the dimensionality from 2NM to (N + 1)M. The most noticeable difference is the larger number of states when ensembles are used (Fig. 2b, inset). Despite the larger number of states, plotted on the same energy scale, a non-diminishing gap between the ground and excited-state maintained for the ensemble case (Fig. 2b, main figure). This occurs due to the larger energy scale of the ensemble Hamiltonian by a factor of N, which at least partially compensates for the larger number of states. Many of these additional states are logically equivalent states in the sense of Eq. (6). For example, we label the states at λ = 1 in terms of the eigenstates \({S}_{i}^{z}\). In the qubit version, the two lowest states have a spin configuration of \(({\sigma }_{1}^{z},{\sigma }_{2}^{z},{\sigma }_{3}^{z})=(-1,-1,-1)\) and (+1,−1,−1), respectively. In the ensemble version with N = 5, of the shown states, the lowest 7 states are all logically equivalent to (−1,−1,−1) in terms of (Eq. 6). Such logically equivalent states provide protection against error since they occur with energies in the vicinity of the ground state, and act as a "buffer” before logical errors are induced47,50.

Fig. 2: Energy spectrum and gap energies of the adiabatic quantum computing Hamiltonian.
figure 2

The spectrum of Eq. (3) with M = 3 with (a) N = 1 and (b) N = 5 for parameters J12 = − 0.5, J13 = 0, J23 = − 1, K1 = 0.5, K2 = 0, K3 = 1. The mean-field approximation for the N = 5 is shown as the dashed lines for the ground and first excited state. c The gap energy for the ensemble qubit numbers is shown. d Scaling of the gap energy with N for two values of λ as shown. The mean-field (MF) approximation is shown as the dotted line in c and d. e The minimum gap vs. N for 60 random instances of problems with Nc = 3. The average minimum gap, as well as the largest (best) and smallest (worst) instances, are shown. f Same as e for problems with Nc > 7.

Full size image

Our aim in the AQC will be to keep the adiabatic evolution in the ground state of the ensemble system. Obtaining the ground state and first excited state for the ensemble system, in general, is a numerically intensive task involving a diagonalization within a Hilbert space of dimension (N + 1)M. To see the behavior for large ensemble sizes, it is desirable to have an approximate method of estimating the gap energy that does not require full diagonalization. Mean-field theory provides an accurate estimate of physical quantities for large spin systems. The ensemble nature of the Hamiltonian allows us to extract energies with increasing accuracy particularly for large N. We use a mean-field ansatz wavefunction of the form

$$\left|{{{\Psi }}}_{\,\text{MF}\,}^{(0)}\right\rangle =\mathop{\prod }\limits_{i=1}^{M}{\left|0,{\theta }_{i}\right\rangle }_{i},$$
(9)

where we define a Fock state of N spins all aligned in the same direction as

$${\left|0,\theta \right\rangle }_{i}=\mathop{\prod }\limits_{n=1}^{N}\left(\cos \frac{\theta }{2}{\left|0\right\rangle }_{ni}+\sin \frac{\theta }{2}{\left|1\right\rangle }_{ni}\right),$$
(10)

which is the maximal positive eigenstate of the rotated spin operator \({\tilde{S}}^{z}=\sin \theta {S}^{x}+\cos \theta {S}^{z}\). We note that a similar mean-field ansatz was used in past works to analyze the N = 1 AQC Hamiltonian ground state38,40,44,62. We apply the mean-field ansatz by performing a self-consistent procedure to obtain the parameters θi. This is equivalent to taking expectation values of the Hamiltonian (Eq. 3) with respect to Eq. (9) and optimizing for θi (see Supplementary Note 2). From the discussion relating to the logically equivalent states, a suitable mean-field ansatz for the first excited state consists of a spin-wave state where one qubit per ensemble is flipped

$$\left|{{{\Psi }}}_{\,\text{MF}\,}^{(1)}\right\rangle =\mathop{\sum }\limits_{k=1}^{M}{\psi }_{k}{\left|1,{\theta }_{k}\right\rangle }_{k}\mathop{\prod}\limits_{i\ne k}{\left|0,{\theta }_{i}\right\rangle }_{i},$$
(11)

where we have defined

$$\left|1,\theta \right\rangle ={\tilde{S}}^{x}{\left|0,\theta \right\rangle }_{i}.$$
(12)

and \({\tilde{S}}^{x}=-{S}^{x} \cos \theta +{S}^{z} \sin \theta\) which creates a spin-flip in the \({\tilde{S}}^{z}\)-basis. To apply the mean-field ansatz (Eq. 11) we diagonalize an effective Hamiltonian in the ψk coefficients and take the lowest energy state (see Supplementary Note 2). We note that the mean-field theory is only expected to work in the regime with N ≥ Nc, since the first excited state is taken to be of the form (Eq. 11), which has a spin configuration that is one spin-flip away from the ground state.

The results are shown in Fig. 2b, where the mean-field estimates (dashed lines) are compared to the exact results. We see that excellent agreement in the energies of the states is obtained for all values of the adiabatic parameter λ. In Fig. 2c we plot the exact gap energy for various N in comparison to the mean-field calculation. Figure 2d shows the convergence of the energies towards the mean-field results with N at various intermediate values of λ. The mean-field results correspond to the limit N → ∞, and the exactly calculated gaps for various N rapidly approach the mean-field result.

The results of Fig. 2a–d were for a particular problem instance. What is more meaningful is to study the performance of the scheme for a variety of different problem instances so that the overall behavior can be assessed. We find that the behavior is rather different depending upon whether N < Nc or N ≥ Nc, due to the different nature of the first excited state. We study the two regimes separately by choosing problem instances where Nc occurs relatively early (Nc = 3, Fig. 2e) or late (Nc > 7, Fig. 2f). In Fig. 2e, we show the average, best, and worst scaling of the minimum gap for problems with Nc = 3, such that most of the N-dependence are in the N ≥ Nc regime. The best and worst scalings are defined as the largest and smallest difference in the gap comparing the qubit and N = 7, the largest ensemble size calculated. We find that the minimum gap increases with N on average. Combined with the logically equivalent buffer states in the region of the ground state, we expect that the AQC performance should improve for these cases.

For the cases with Nc > 7 (the N < Nc regime), we see more mixed results (Fig. 2f). The average scaling tends to still improve with N, but there are some cases where the minimum gap becomes significantly worse with N. In such cases we expect that the AQC performs poorly. We note that the small values of N considered here are due to limitations in our numerical simulations. We thus expect that realistic ensemble sizes would satisfy N≥Nc, where the scaling is more favorable.

This may, at first glance, seem to be a counter-intuitive result, since one might expect that with a larger N the system should behave more classically. However, it can be seen that in both Eqs. (4) and (5) the magnitude of the elementary excitation does not diminish as N grows, since it is always a discrete Hamiltonian. Thus, the gap does not diminish even for N → ∞, and AQC can be performed with macroscopically sized ensembles.

Performance with adiabatic evolution

We now directly time-evolve the AQC Hamiltonian and demonstrate its performance. We use a linear annealing schedule λ(t) = t/τ and examine the final occupation probability of the states at time t = τ. First, examining the case without decoherence, we vary the sweep time τ for the same problem instance as shown in Fig. 2b. In Fig. 3a, we see that as expected for sufficiently long τ, most of the population is concentrated in the ground state. Diabatic excitations are seen for smaller τ depleting the population in the ground state. Here we also indicate the distribution of the energy levels that are logically equivalent to the ground state (green levels) and error states which have a different configuration (red levels) according to the majority vote encoding (Eq. 6). We see that the diabatic excitations tend to distribute the probability with a tendency to excite the lower energy states. Since the logically equivalent states are also in the low-energy range, this shows that an effective buffer is provided by the encoding, protecting the computation as long as the diabatic excitations are within the logically equivalent states.

Fig. 3: Occupation probability of the final state after time-evolving the adiabatic quantum computation.
figure 3

The occupation probability is plotted as a function of the (a) sweep time τ with N = 5 and Γ = 0; (b) Sz decoherence rate Γz with N = 5; (c) Sx decoherence rate Γx with N = 5; (d) Sz decoherence rate Γz with N = 1. The same parameters as Fig. 2b are used M = 3J12 = − 0.5, J13 = 0, J23 = − 1, K1 = 0.5, K2 = 0, K3 = 1, and τ = 100. The probability distribution is normalized such that the maximal probability takes a value of 1 for each τ or Γx,z, for the clarity of the figure. For each plot, the spectrum of energy levels is labeled to the right of the plot. Green levels indicate logically equivalent states to the ground state, red levels indicate logical error states.

Full size image

The addition of decoherence has a qualitatively similar effect to diabatic excitations, as can be seen from Figs. 3b–d). We numerically evolve a master equation in the presence of Markovian Lindblad dephasing in the Sz and Sx basis and obtain the final probability distributions for various decoherence rates (see “Methods” section). Both the Sz-dephasing and Sx-dephasing is found to have a similar effect, with energy levels in the lower region being populated for stronger decoherence rates (Figs. 3b, c). Again, due to the fact that the logically equivalent states are dominated towards the lower end of the energy spectrum, this shows that the effect of the decoherence will also be to initially excite the logically equivalent states.

It is interesting to compare the distribution to the unencoded bare AQC Hamiltonian N = 1 (Fig. 3d) to the N = 5 encoded case (Fig. 3b), plotted on the same energy scale. We see that for small decoherence rates, the effect of the dephasing is similar in terms of the width in the energy of the probability distribution. The benefit of the encoding scheme in the low decoherence region is evident from the fact that the first logical error state has much higher energy in the N = 5 case rather than the N = 1 case (9 and 3 units above the ground state, respectively). Summing the probability of all the states that are logically equivalent to the ground state (green levels), the ensemble case has a higher success probability. However, for larger dephasing rates, the distribution for the ensemble case becomes broader, and the error suppression is less effective. We, therefore, expect that as long as the decoherence rate is below a threshold, the overall logical errors can be effectively suppressed. We note that while we illustrate our results with a single problem instance, we have examined a variety of cases and seen consistently similar results for other problem instances.

We now examine the dependence of the logical errors at the end of the AQC with the ensemble size N. In this case, it is illustrative to examine another problem instance, corresponding to the ferromagnetic Hamiltonian Jij = − 1(1 − δij) with a bias field Ki = K. For N = 1 with K > 0, the ground state is σi = − 1 and the first excited state is σi = + 1 for all i. The energy landscape corresponds to one global minimum and one local minimum separated by a potential barrier consisting of all the remaining states. For the ensemble case and in the regime N≥Nc, the ground state has the same logical configuration \({S}_{i}^{z}=-N\) but the first excited state is a single spin-flip of the configuration \({S}_{k}^{z}=-N+2,{S}_{i}^{z}=-N\) ∀ i ≠ k (see Supplementary Note 3). The error probability is defined as 1 − (success probability), where the success probability is defined to be the total probability of all states that are logically equivalent to the ground state.

 Figure 4a, b shows the error probabilities for Sz and Sx dephasing, respectively, for the ferromagnetic instance with different values of K. We observe that the error probabilities increase initially, but start to strongly decrease as N is increased further. While lower error probabilities are expected from a larger gap with larger N as observed in Fig. 2e, the decrease in error is much more than would be expected from this. The suppression of the error with increasing N attributed to the fact that above Nc, the first excited state becomes a logically equivalent state, which is one of the characteristics of the N≥Nc regime (Fig. 1a). As the effect of the dephasing is to excite low-lying energy states, the excitation of the first excited state no longer becomes a logical error, suppressing the total error. Such an improvement is consistent with the analysis in ref. 58 where the encoding was shown to give effectively give lower temperatures, and ref. 59 where phase transitions in the model were shown to be weakened.

Fig. 4: Error probabilities in the time-evolved adiabatic quantum computation.
figure 4

a, b Error vs. N for the ferromagnetic model with M = 3, Jij = − 1(1 − δij) and Ki = K, and τ = 100. The critical ensemble size for each parameter is: Nc(K = 0.1) = 14, Nc(K = 0.2) = 8, Nc(K = 0.3) = 5. Dephasing in the basis (a) Sz-basis with Γz = 10−4 and (b) Sx-basis with Γx = 10−4 are calculated. c, d Averaged error vs. N for various τ and M = 3 for (c) 60 problem instances with Nc = 3; (d) 60 problem instances with Nc > 7. Averaged error vs. N for various τ and M = 3 including dephasing of rate Γ = 10−4 for the same problem instances with (e) Nc = 3; (f) Nc > 7.

Full size image

The phenomenology of an error increase in the region for small N but decreasing error for large N is also observed in other problem instances. We randomly generate problem instances and calculate the error probability of the ensemble encoding as a function of various ensemble sizes N and sweep times τ in the absence and presence of decoherence. For each generated problem instance, we calculate the critical Nc, and again group the instances according to whether Nc occurs relatively early (specifically Nc = 3) or late (Nc > 7). Examining each case separately, we can study the performance of the ensemble scheme in the N ≥ Nc or N < Nc regime, respectively. Figure 4c–f shows the error probability of finding the state of the system at the end of adiabatic evolution, averaged over a set of fixed 60 randomly chosen problem instances with Nc in each range.

First examining problems with Nc = 3 (the N ≥ Nc regime) without decoherence, we observe that for sufficiently long τ the error probability strongly decreases with N (Fig. 4c). We again attribute this to the presence of logically equivalent excited states in the vicinity of the ground state. Although we only simulated relatively small system sizes due to numerical limitations, we expect that for larger N the trend will continue towards lower errors as the gap energy approaches the mean-field value corresponding to N → ∞. For short sweep times τ, the errors also improve with N, although with a smaller gradient. We attribute this behavior due to the sweep times being in a diabatic regime, such that the system is not maintained in the vicinity of the ground state, which involves high energy excitations. We note that there did exist rare examples where the error probability scaled badly with N, due to the particular structure of the energy spectrum. However, the occurrence of these poorly scaling examples was so rare that they did not impact the average to a significant extent.

For problems with Nc > 7 (the N < Nc regime) and no decoherence, the error tends to increase with N (Fig. 4d), despite the fact that the average minimum gap increases with N, as seen in Fig. 2f. We have examined the individual cases and confirmed that for particular cases where the gap increases with N, the error decreases with N as expected. The reason that the average error increases are that the cases with poor gap scaling in Fig. 2f tend to have nearly zero success probability, which reduces the average, considerably. Thus, the results for problems in the N < Nc regime are mixed and depend very much upon the particular problem instance of whether the gap increases or decreases. We do, however, note that these problem instances themselves are rather rare, as seen in Fig. 1c, in comparison to the more common N > Nc case. When averaged over all random problem instances, errors tend to decrease with N.

Calculations adding Sz-dephasing to the AQC scheme are shown in Fig. 4e, f. We find generally the same behavior of the error with N when decoherence is introduced, but with a higher error probability overall, as expected. For the Nc = 3 case shown in Fig. 4e, we see that there is a similar improvement of the error with N as the zero decoherence case. The new feature here is that there is an optimum sweep time beyond which the error probability starts increasing again. This can be simply explained by noting that the AQC must be performed within the decoherence time available to the computation. Beyond the optimal time, the performance starts to degrade, therefore there is a trade-off between maintaining adiabaticity and working within the decoherence time63,64. For the Nc > 7 cases in Fig. 4f, we again see the error increase with N, which is attributed again to particular instances where the minimum gap decreases with N.

We note that in Figs. 3 and 4 we have used collective dephasing with respect to the Sz and Sx operators, which is one of the decoherence channels for the atomic ensemble implementation65,66,67,68. For other implementations, individual qubit dephasing may be a more relevant decoherence model. We generally expect similar results with individual qubit dephasing, since elementary perturbations due to decoherence produce similar transitions in the energy spectrum of the Hamiltonian. Our calculations based on the ferromagnetic model confirm that qualitatively similar results are obtained with an individual qubit dephasing model. The primary difference, in this case, is the addition of non-symmetric states, which can also act as buffer states to protect the ground state (see Supplementary Note 4).

We finally comment on the presence of entanglement during the AQC sweep. The mean-field wave function as given in Eq. (9) takes the form of a product state of spin coherent states on each ensemble. This may suggest that there is no entanglement between the ensembles during the adiabatic evolution. In fact, entanglement is typically present during the evolution, due to the \({S}_{i}^{z}{S}_{j}^{z}\) interaction in the AQC Hamiltonian. The mean-field ground state is merely an approximation to the true ground state, which in fact typically contains entanglement. We have explicitly calculated entanglement for small ensemble sizes (see Supplementary Note 5). The presence of entanglement is consistent with past works studying the robustness of entanglement in the presence of decoherence60,61,66,69,70,71,72,73. The factor of 1/N multiplying the \({S}_{i}^{z}{S}_{j}^{z}\) terms makes the type of entangled state of a robust type as discussed in ref. 66. We, therefore, expect that the entanglement should survive for macroscopic ensembles within the decoherence window. This can be contrasted to other ensemble-based approaches such as in liquid-state NMR74,75, where the entanglement collapses to zero.

Experimental implementation

We briefly describe a potential experimental implementation that realizes the Hamiltonians (Eq. 4) and (Eq. 5). Neutral atom ensembles, consisting of either thermal atomic ensembles or cold atoms, are a strong candidate for realizing a single spin ensemble \({S}_{i}^{x,z}\) and 60,68,70,71,76. The individual qubits of spin \({\sigma }_{i,n}^{x,z}\) within the ensemble are realized by the internal hyperfine ground states of the atoms. For example, for 87Rb the ground states F = 1, mF = − 1 and F = 2, mF = 1 are clock states such as their energy separation is insensitive to magnetic field fluctuations68. Thermal atomic ensembles are trapped in paraffin-coated glass cells which prolong the coherence time of the internal states70. In this case, each glass cell acts as one of the ensembles spins \({S}_{i}^{x,z}\), such that M glass cells are prepared. Such a multi-ensemble system was realized in ref. 77 where 225 locally addressable atomic ensembles were created, as well as entanglement generation between 25 ensembles78. Another approach is to have a multi-trap atom chip system, such as that proposed in ref. 65,76. Atom chips are a flexible platform for producing magnetic traps, where atoms can be cooled to quantum degeneracy. One advantage of cooling is that long coherence times can be achieved (up to 1 min in ref. 79), for atomic gases just above the critical temperature for Bose-Einstein condensation.

To realize Hamiltonians (Eq. 4) and (Eq. 5), we require both ensemble-ensemble effective interactions and single ensemble control. This can be produced by using optically mediated methods, where off-resonant lasers produce an effective \({S}_{i}^{z}{S}_{j}^{z}\) type interaction69,80,81. An alternative for cold atom systems is to produce the interactions by taking advantage of the non-linear interactions between the atoms using state-dependent forces82. Optically mediated methods are flexible as they are able to produce remote entanglement, whereas the interaction methods require bringing the ensembles into close proximity. Adjusting the interaction parameters allows one to realize a controllable Jij matrix between the ensembles. The Hamiltonian (Eq. 5) can be achieved by microwave/radio frequency transitions. In the case of 87Rb, this is realized by a two-photon transition between the clock states67. Unlike optical frequencies which can be focused on a single ensemble, microwave/radio frequency transitions are less easily directed and are applied on all ensembles in the same vicinity, which is sufficient for the \({S}_{i}^{x}\) Hamiltonian. Finally, the \({S}_{i}^{z}\) terms (Eq. 4) could be realized by a state-dependent potential83 or an ac Stark shift to optically shift the energy levels76. We note that other physical systems could potentially implement the Hamiltonians (Eq. 4) and (Eq. 5), for example, ensembles of NV-centers84,85,86 are another possibility.

One of the attractive features of the AQC Hamiltonian (Eq. 4) and (Eq. 5) is that it is written completely in terms of the total spin of the ensembles of qubits. This means that only global control of the ensembles—involving the collective spin operators \({S}_{i}^{x,z}\) —rather than the control of individual qubits \({\sigma }_{i,n}^{x,z}\) is necessary. Alleviating the need for individual qubit control in realizing an error-protected quantum computer is a significant simplification since most quantum error-correction schemes require rather sophisticated degrees of quantum control87,88. The resources required to scale up the degree N of the repetition code would be considerably less in an approach involving only collective spin control. For example, in atomic ensembles on atom chips, one typically deals with atom clouds of the order N ~ 103 and 61,67, and in paraffin coated glass cells the number of atoms is N ~ 1012 and 60, which are all manipulated in parallel within the experiment. In systems where each of the copies must be implemented with individual qubit control, the resources are consequently greater. In comparison, past demonstrations on the D-wave machines have been in the region of N ~ 1050,56,57. While related approaches have been considered in the context of all-to-all encodings such as nested QAC56,57,58, the chimera graph hardware that was used necessitates a further minor embedding step which is not required in our case.

An important issue when using macroscopic ensembles is sensitivity to decoherence. One might naively expect that atomic ensembles containing up to N ~ 1012 atoms are very sensitive to decoherence. In fact, the sensitivity of the quantum state to decoherence is a highly state-dependent process60. For example, cat states are extremely sensitive to decoherence, but spin coherent states are more robust66,89. Hence the important consideration is the type of states that are generated, and one should be careful that they are robust in the presence of decoherence. It was shown in ref. 66 that the \({S}_{i}^{z}{S}_{j}^{z}\) interaction produces states which are robust in the presence of decoherence as long as the interaction timescale is of order 1/N. In this respect, the factor of 1/N in front of the two-ensemble interaction (Eq. 4) is beneficial as it implies highly decoherence-sensitive states are not generated. Another indication of this is that in the N≥Nc regime, the mean-field state (Eq. 9) is a tensor product of spin coherent states, which are known to be robust in the presence of decoherence. This explains the good performance of our proposed scheme even in the presence of dephasing, as shown in Fig. 4.

Further details of the implementation with neutral atom ensembles can be found in refs. 65,76,90 and Supplementary Note 6.

Discussion

In this study, we investigated a formulation of AQC where qubit ensembles are used instead of qubits, and the ensemble Hamiltonians (Eq. 4) and (Eq. 5) are adiabatically evolved. We have found that finding the ground state of (Eq. 4) is an equivalent problem to the original qubit problem Hamiltonian (Eq. 1). The main difference between the ensemble and qubit problem Hamiltonians is that the ensemble version introduces many logically equivalent states as defined in (Eq. 6) with similar energies to the ground state. The introduction of these states is beneficial for AQC since the occupation of these states does not cause a logical error, and provides a buffer against diabatic excitation. We found that there are two important regimes with respect to N, depending on the character of the first excited state, summarized in Fig. 1a. In the regime with N ≥ Nc, we find that the minimum gap energy increases, and the ground state is logically protected, leading to a reduced error probability in the AQC. In the regime with N < Nc, we obtain mixed results, where despite the average minimum gap increasing, the AQC scales on average poorly. This was due to the particularly poor performance of the cases where the gap decreases and can be attributed to the lack of logical protection of the ground state. For large ensemble sizes such as that realized with atomic ensembles, all but a minority of problems should satisfy N ≥ Nc, where the ground state is logically protected.

We thus find that AQC with ensembles should perform well in a great majority of cases for large N. One may find it surprising that it is possible to perform AQC at all with ensembles of qubits, even in the limit of N → ∞. The first key point that allows for the ensemble version to still work is that the discrete nature of the Hamiltonian is preserved under (Eq. 4) and (Eq. 5). Thus, although the energy of the full space can be viewed as being quasi-continuous as shown in Fig. 1a, this is only because space is being viewed in rescaled variables \({x}_{i}={S}_{i}^{z}/N\). Physically, the magnitude of the spins is also growing with N, which preserves the energy gap. From a resource point of view, one may argue that many more physical qubits are being used. However, we take the point of view that the relevant resource is the complexity of the experiment control when dealing with an ensemble as compared to a single qubit. For many implementations, the effort required for controlling an ensemble is no more than that of a single qubit. For instance, if performing a single qubit operation on an atom is performed by a laser pulse, then the equivalent ensemble operation is to apply the pulse to the whole ensemble. This is typically not an operation that costs N times more since one can illuminate the whole ensemble with the same laser, i.e., it is parallelizable. Thus, as long as the operations for the qubit operators \({\sigma }_{i}^{x,y,z}\) can be performed with a similar experimental overhead to ensemble operators \({S}_{i}^{x,y,z}\), then implementing the ensemble and qubit version of the AQC Hamiltonians should be of comparable complexity.

One may also be concerned that the use of ensembles may be problematic since they could be extremely sensitive to decoherence, owing to their macroscopic nature. The sensitivity of qubit ensemble states has already been investigated in numerous works, see for example refs. 66,90,91,92. The main point here is that the fragility of the quantum states is state-dependent: while Schrodinger cat states are extremely sensitive, spin coherent states are generally quite robust. This is what has allowed the experimental realization of macroscopic quantum states, such as those performed by Polzik and co-workers60,70,71. The form of the mean-field ground and excited states suggest that the ensemble version of AQC can also be robust for the same reasons. The ground state (Eq. 9) is nothing but a set of spin coherent states, and (Eq. 11) is a spin-wave excitation on the ground state. Spin-wave states are also relatively robust and have been already demonstrated experimentally93. Therefore, as long as the ensemble size is such that N ≥ Nc, we believe that it is reasonable to expect that the scheme works even in the presence of decoherence. If it is the case that N < Nc, it is less clear what the decoherence properties are since the nature of the state is not yet understood. We nevertheless observe that in some cases the minimum gap can increase, making the ensemble framework viable. While we have not been able to exactly characterize the cases that are most susceptible, we also have not seen any correlation with classically hard instances of combinatorial problems (see Supplementary Note 4). Considering that these are a small fraction of the full problem set for large N, we find that in most cases the ensemble framework successfully performs error suppression via the duplication of the quantum information. This is consistent with other approaches using repetition codes with AQC50,51,52,53,54,55,56,57,58,59.

Another direction that could be further investigated is the use of energy penalty terms, which have been shown to be beneficial in several works50,56,57,58,59. This would involve the addition of terms to the Hamiltonian of the form \({({S}_{i}^{z})}^{2}\) to induce a ferromagnetic interaction within the qubits in the ensembles. This is expected to further improve the performance, where there is an optimal strength of the interaction. Producing such an interaction in the case of atomic ensembles has been experimentally and theoretically investigated67,69, and is compatible with the general framework of our approach since it is based on collective operations of the ensemble. We will leave these topics for further investigation as future work.

Methods

Numerical simulation

To examine the performance of the ensemble version of AQC, we performed both a pure state evolution and mixed state evolution of Hamiltonian (Eq. 3). We use a linear annealing schedule λ(t) = t/τ and examine the final occupation probability of the states at time t = τ. For the case that we include decoherence, we consider Markovian Sz-dephasing and Sx-dephasing. This is particularly relevant for implementation with atomic ensembles, where the coupling to the ensemble spins occur in a collective manner65,66,67,68. We use the master equation94

$$\begin{array}{l}\frac{d\rho }{dt}=i[\rho ,H]-\frac{{{{\Gamma }}}_{z}}{2}\mathop{\sum }\limits_{n=1}^{M}[\rho {({S}_{n}^{z})}^{2}-2{S}_{n}^{z}\rho {S}_{n}^{z}+{({S}_{n}^{z})}^{2}\rho ]\\ \qquad-\frac{{{{\Gamma }}}_{x}}{2}\mathop{\sum }\limits_{n=1}^{M}[\rho {({S}_{n}^{x})}^{2}-2{S}_{n}^{x}\rho {S}_{n}^{x}+{({S}_{n}^{x})}^{2}\rho ],\end{array}$$
(13)

where Γz,x is dephasing rates and H is the Hamiltonian (Eq. 3). Starting from the eigenstate of the initial Hamiltonian we solve the master equation numerically for the combined adiabatic and dephasing evolution. The performance of the AQC is then evaluated through the probability of finding the state of the system in the ground state at the end of the adiabatic evolution. Qutip was used for the simulations95.