We consider a chain of ions confined in a linear Paul trap. The spin-1/2 Ising variables appearing in H_final can be encoded by restricting the internal atomic dynamics to only two hyperfine sublevels. Single-qubit rotations allow for a preparation with high fidelity of any desired initial product state [13], including the completely polarized ground state of H_init [6–8, 10–13]. Similarly, the projection of the spin-1/2 variable on any coordinate direction can be measured to high precision by appropriate single-qubit rotations followed by stimulated fluorescence measurements, allowing the reconstruction of arbitrary correlation functions [7, 13]. The synthetic magnetic fields h^z_i and h^x appearing in H_final and H_init can also be simulated in a straightforward fashion via AC-Stark shifts or by resonantly driving the transition between the two sublevels [40].

Ising interactions can be simulated by coupling the ions via laser [40–42] or microwave radiation [43] to the collective vibrational modes. In the limit of off-resonant radiation, one can disentangle the collective vibrations and the atomic Ising variables via a canonical transformation and obtain an effective Hamiltonian for the Ising variables alone. For example, when employing a Mølmer–Sørensen scheme [42], the Ising spins are subject to the effective interactions (see Appendix A) (4)Jij=−ℏ|Ωi||Ωj|4∑qηiqηjqδq, where η_iq is Lamb–Dicke parameter including the amplitude of vibrational mode q at ion i, δ_q is the detuning of the laser relative to the vibrational mode q, and |Ω_i| is the absolute value of the laser's Rabi frequency at ion i. For the derivation of Equation (4), we assumed the off-resonance condition η_iq|Ω_i| ≪ δ_q. The addition of a transverse field, required for our choice of H_init, induces additional couplings to the phonon modes in higher order, which, however, can be neglected as long as h^ν_i ≪ ħδ_q, ν = x, z. The precise form of the interactions depends on the coupling scheme used, but the qualitative behavior is similar for other schemes [40, 43].

As shown, e.g., in Britton et al. [10], Islam et al. [8], Jurcevic et al. [13], when transmitted via the transverse phonon modes the Ising interactions follow an approximate power law as a function of the distance, J_ij ~ 1/|i−j|^α, with a decay coefficient 0 ≤ α ≤ 3 that can be adjusted using the detuning δ_q (where larger detuning entails larger α). When employing the longitudinal phonon modes, the large spacing of vibrational frequencies restricts the range to α ≈ 0. By expressing the local Rabi frequencies as multiplies of a reference frequency |Ω|, i.e., |Ω_i| = F_i|Ω|, we can thus write (5)Jij=JFi Fj|i−j|α, with J setting the overall energy scale. Note that the power-law decay is exact only in the limits α = 0 and α = 3. As shown recently in Karp [44], the interactions in the intermediate regimes are better described by an exponential tail added to a dipolar decay. Additional modifications may arise from an uneven spacing of ions as is usual in linear Paul traps or from inhomogeneous laser profiles. For small chains, a general power-law behavior is, however, a good approximation of the system dynamics [13].

Until now, we considered a coupling to a single set of vibrational modes. However, richer models can be quantum simulated by employing the fact that in a linear Paul trap the three orthogonal directions of vibration decouple [40]. In this case, by using different lasers to couple to each orthogonal set of modes, it is possible to generate three independent Ising interaction terms, which we label ℓ = 1, 2, 3. Additionally, one may induce synthetic fields h^z_i in an ion-dependent fashion by, e.g., resonantly driving the qubit transitions with beams derived from an acousto-optic deflector [13]. It is thus possible to realize in trapped ions—using available technology—the following Ising Hamiltonian:

Experiments thus far—aiming at mimicking translationally invariant many-body models—made great efforts to engineer systems as homogeneous as possible, i.e., F^(ℓ)_i = 1 and h^z_i = h^z, ∀i [7, 8, 10, 12, 13]. For the same reason, few theoretical works have considered a site-dependent tunability of couplings [38, 45–47] or fields [48]. In contrast, our purpose of engineering NP-complete problems in Hamiltonian Equation (1) requires a certain degree of programmability of the interaction matrix elements J_ij and fields h^z_i (see e.g., Lucas [16]). In Table 1, we present several NP-complete problems that can be implemented with interactions of the form given in Equations (5) and (6), i.e., without requiring a full programmability of interactions. These problems thus become amenable to current trapped-ion technology by adjusting the local laser intensities |Ω_i|².

The NP-complete problems implementable in this way include the number-partitioning problem, the integer knapsack problem [16], and the Coulomb glass problem, which—as shown in Anjos et al. [49]—can be mapped to a class of opti-cut problems. We summarize the description of the necessary Hamiltonian parameters in Table 1. In the numerical studies presented in Section 5, we focus on the example of the Coulomb glass problem, since it is closest to current experiments [12, 13]. The Coulomb glass [51] contains a random magnetic field h^z_i ∈ [−ϵ_z, ϵ_z] and Ising interactions that are homogeneous and long-ranged, α⁽¹⁾ = 1 and F⁽¹⁾_i ≡ 1. Usually, one imposes the constraint of vanishing total magnetization, which can be realized by mean-field-like interactions (α⁽²⁾ = 0) with F⁽²⁾_i ≡ F⁽²⁾ ≫ 1. Defining the strength of the constraint as V ≡ J(F⁽²⁾)² ≫ J, the total Hamiltonian is then (7)Hfinal=J∑i≠jσizσjz|i−j|+∑ihizσiz+V∑i≠jσizσjz, with J, V > 0. The large frustration of the Ising interactions counteracts ordering tendencies and, for ϵ_z = (J), leads to a strong competition between the interactions and the random local fields. This competition signs responsible for its computational complexity. Currently, it is not clear whether the Coulomb glass exhibits a spin-glass transition at nonzero temperature. This, however, does not necessarily imply a reduced complexity of the ground state of the model. Certain random Ising models on Chimera graphs, such as implemented in the D-Wave device, for example, do not show spin-glass transitions at nonzero temperature, but finding the ground state of the model is still NP-hard. Additionally, it may be that a non-perfect power-law decay, for example the more realistic dipolar plus exponential shape [44] affects the complexity of the model.

A fundamental question in AQO is how the optimization of problems such as given in Table 1 performs in the presence of decoherence. This question is of particular importance for current implementations of AQO devices such as the D-Wave machine where decoherence rates are large compared to the annealing times. Current trapped-ion setups are able to suppress noise to very low values, currently with time scales around 10/J [13]. This opens the possibility to systematically study the performance of AQO under varying degrees of purposefully engineered decoherence. A scheme to engineer different types of classical noise in a controlled fashion has recently been discussed and implemented in Soare et al. [52].

Here, we describe how classical white noise can be engineered by a conceptually simple extension of the steps leading up to Equation (4). Detailed derivations can be found in the Appendix C. Here, we sketch the main idea, which is based on the addition of randomly fluctuating synthetic magnetic fields. We assume that the rate of these fluctuations is much larger than any other parameter in the Hamiltonian, so that we can treat the annealing parameter A as quasi-constant. The total Hamiltonian then takes the form (8)H=AHinit+(1−A)Hfinal+∑i[ξix(t)σix+ξiz(t)σiz] , where we choose the noise in two spin directions ν = x, z, and we take ξ^ν_i(t), to be independent Ornstein–Uhlenbeck processes, 〈〈ξiν(t)ξjμ(0)〉〉=ℏΓδν,μδi,jℏbe−b|t|. Here, 〈〈·〉〉 denotes the average over noise realizations, and for simplicity we take the rate Γ and the bandwidth b to be independent of ν = x, z. The fluctuating fields can be realized using the same experimental techniques as for their quasi-static counterparts described at the beginning of Section 2.1. The additional requirement of fast fluctuations that are uncorrelated between different sites could be solved via a simultaneous addressing of the individual ions using several laser beams derived from an acousto-optic deflector [13].

If the bandwidth b of the fluctuations is much larger than the internal energy scales of the ions and the coupling to the phonon modes, b ≫ h^ν_i/ħ, η_iq|Ω_i|, we can average over the noise and obtain a master equation describing classical white-noise dephasing, (9)∂ρ∂t=−iℏ[H,ρ]+Γ∑iD[σi(x)]+Γ∑iD[σi(z)] , where D[X] = 2XρX^† − ρX^†X − X^†Xρ is the Lindblad super-operator. In the derivation of Equation 9, we neglected additional couplings to the phonon modes that appear in higher order and can be dropped provided b ≫ η_iq|Ω_i|,h^ν_i/ħ. These conditions complement the requirements for the derivation of AH_init + (1 − A)H_final, which read δ_q ≫ η_iq|Ω_i|, h^ν_i/ħ. In the Appendix C, we discuss that they can all be fulfilled for realistic experimental parameters.

Hence, adding fluctuating fields with bandwidth much larger than other relevant energy scales allows to systematically engineer white noise in trapped ions. By inducing such fluctuations on the longitudinal or transversal fields, it is possible to study different types of Lindblad super-operators.

A convenient, well-established measure for entanglement of closed systems is the entanglement entropy, which is defined between two subsystems A and B. The entanglement entropy is the von Neumann entropy of the reduced density matrix ρ_A = tr_Bρ of A obtained by tracing out the complement B,

In this work, we study the maximum of the entanglement entropy considering all possible pairs of particles as subsystem A and the remaining N − 2 particles as subsystem B. This pair-block entanglement entropy maximized over all pairs is denoted as S^max₂(t). The largest pair-block entanglement entropy reached during the sweep is S^max₂ = max_t [S^max₂(t)].

A convenient mixed-state entanglement measure is the logarithmic negativity between two parts C and D of a larger subsystem A (A = C ∪ D). It is defined as [53, 54] where ||X||=trX†X is the trace norm. ρTCA is the partial transpose of the reduced density matrix ρA with respect to subpart C, given by 〈dc|ρTCA|d′c′〉 = 〈dc′|ρA|d′c〉 where the c's and d's form a basis for C and D, respectively. In the following, we consider the logarithmic negativity between subsystems C and D consisting of a single site each, and maximize again over all combinations.

While entropy and negativity are measures for the entanglement between two subblocks of the system, we are also interested in multi-particle entanglement. This can be measured with the Fisher information F_Q, which not only quantifies entanglement, but also provides a lower bound on the number of particles that are entangled [55, 56]: if a state is k-producible—meaning that it can be written as a product of states involving less than or equal to k particles each—the Fisher information satisfies

Here, N is the total number of spins and ⌊ · ⌋ denotes the floor operation. If F_Q > N, this inequality implies entanglement between an increasing number of particles with increasing F_Q, where F_Q close to the maximal value of N² corresponds to N-particle entanglement. In quantum annealing devices, the multi-particle nature of entanglement may be a crucial aspect, since in a spin glass, for example, the lowest-energy configurations are typically vastly different. Thus, many spin flips should be necessary to tunnel from one local energy minimum to the next. If one assumes that all of these spins have to be coherent during the tunneling process, multi-particle entanglement may become the critical phenomenon for the efficiency of AQO protocols.

While the above entanglement observables are powerful tools to theoretically understand the nature of entanglement, they require elaborate measurements in experiment. Therefore, it is desirable to find entanglement witnesses that can be deduced from simple observables.

Here, we introduce an optimal local entanglement witness W that allows one to measure entanglement in closed as well as open systems. It can be obtained from two-point spin correlations, (13)W=max{0,−minf→[Wf→] }, with where f_i ∈ ℂ is subject to the normalization ∑_i|f_i|² = 1. This constitutes a generalization of the witnesses introduced in Krammer et al. [57], Cramer et al. [58], De Chiara and Sanpera [59], Hauke et al. [60]. Going beyond Hauke et al. [60], we analyze here the “optimal” witness, obtained by minimizing W_f→ over all f→. In practical terms, this minimization is extremely simple: by introducing a Lagrange multiplier for the condition of normalization ∑_i|f_i|² = 1, it becomes equivalent to finding the smallest eigenvalue of the correlation matrix _ij = ∑ν〈σ^iνσ^jν〉. This procedure can be done via post-processing of the measured two-point correlation functions, which are simple observables in trapped-ion experiments [6–13]. Thus, it is straightforward to determine the optimized entanglement witness experimentally. Since no cross-correlations of the form 〈σ^ixσ^jy〉 are required, it needs only three series of measurements for the x, y, and z correlations respectively, and is therefore easily scalable.

The above entanglement witness gives a lower bound on the best separable approximation (BSA) [58], an entanglement measure that is applicable also to mixed states. The ease of measuring this witness, both in theory and experiments, may make it an important entanglement quantifier, especially for large systems.

For computing the dynamics in closed systems, we employ exact diagonalization with Krylov time evolution [63] for N = 16 spins, and with about 100 disorder realizations of h^z_i for each value of t_a. Note, that the focus of this work is the entanglement during the optimization with restricts the number of spins to relatively small sizes.

Let us first illustrate the dynamics of the entanglement in typical sweeps. Figure 2 depicts the entanglement dynamics and success probability of a few instances of randomly chosen magnetic fields h^z_i. In these examples, the annealing time is fixed to a moderate value t_a = 500/J, such that for most instances the sweep is nearly adiabatic (here, as in what follows, we set ħ to unity). Starting from the separable initial state, entanglement gradually increases up to a maximum around t/t_a ≈ 0.6 − 0.7, i.e., in the vicinity of the equilibrium magnetic phase transition of the homogeneous system. During the final stages of the protocol, the quantum fluctuations introduced via the transverse fields h^x ∑_i σ^x_i become weaker; the entanglement decays while simultaneously the success probability increases. There exist, however, particular “hard” instances where the system is not capable to disentangle, and thus ends up in a superposition of various eigenstates with a lower success probability. The optimal local entanglement witness, which is not an entanglement measure, shows stronger oscillations than the other observables, which results in a peculiar additional peak in the hard instance of Figure 2.

In an ideal, closed system, the only mechanism to lower the success probability is non-adiabaticity due to finite sweeping times. However, in many current implementations of AQO, the efficiency is also reduced by decoherence in the system. In particular, in the D-wave device, decoherence is orders of magnitude faster than annealing times. Trapped ions, on the other hand, operate in a regime where, as we have shown in Section 2, decoherence can be systematically tuned. This ability allows one to study the complex interplay between non-adiabaticity and decoherence in AQO. In this section, we numerically analyze entanglement witness W, negativity E and the Fisher information F_Q in the AQO protocol. We present numerical data from a relatively small system with N = 6 sites where considerable statistics is achievable also for numerically exact open-system calculations. In current experiments, the number of ions can go up to N = 20 [5], which represents an example of a quantum simulator reaching system sizes that cannot be sampled numerically.

We compare various annealing times t_a, ranging from t_a = 4/J to t_a = 256/J, and decoherence rates ranging from Γ = 0 to Γ = 0.01J, depicted in Figure 6. For fast sweeps (t_a = 4/J), the results are independent of Γ, as can be expected since the considered rates are too low to have any effect in this case. In this regime, the final state remains close to the initial equal superposition of all eigenstates of the final Hamiltonian, and the success probability approaches the limit of P = 1/2^N, independent of the rate of decoherence. The entanglement in the final state is then rather uncorrelated from the success probability, which is always low.

In the opposite extreme of slow sweeps (t_a = 256/J), the success probability P is governed almost solely by the decoherence rate, and, contrary to the case of closed systems, has little relation to the final-state entanglement. In fact, in this case of extremely slow sweeps, only at small decoherence does any entanglement as measured through E and F_Q survive, while the witness W vanishes for all studied instances. In the regime of intermediate annealing times, for small decoherence rates the success probability and entanglement are related in a similar fashion as in the fully coherent case, but a large decoherence suppresses both entanglement and success probability.

The optimal witness is rather sensitive to decoherence and does not allow one to obtain detailed information about the process. In contrast, the Fisher information shows rich features and offers additional insight into the nature of the entanglement of the final state. The cases of slow sweeps (t_a = 256/J) and small decoherence rates are particularly interesting for the AQO protocol, since these allow for a large success probability. In these cases, we find that the Fisher information approaches F_Q ≈ N, which indicates entanglement with k ≈ 1, i.e., between individual particles. For faster annealing, F_Q can increase slightly, but still remains much below F_Q = N². We find a similar behavior for the maximum of F_Q assumed during the AQO sweep. This means entanglement between large numbers of particles as measured through F_Q is rare for the studied system sizes.

State-of-the-art trapped-ion experiments have intrinsic decoherence rates on the order of Γ ~ 0.1J, but predominantly only in one quadrature [13]. We have performed numerical tests also for dephasing only and found an improvement of a factor of about 2 compared to the results presented in Figure 6, where we considered the worst case of equal noise in σ^x and σ^z. Thus, for a systematic study of the influence of noise on AQO in such setups, further improvements in terms of decoherence still seem necessary. Also, as is usual with AQO protocols, larger problem instances will show smaller gaps and thus require larger annealing times to achieve similar degrees of success probability. On the other hand, non-linear ramps may improve the performance of the AQO protocol [19, 62].

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.