Improving quantum annealing by engineering the coupling to the environment

A large class of optimisation problems can be mapped to the Ising model where all details are encoded in the coupling of spins. The task of the original mathematical optimisation is then equivalent to finding the ground state of the corresponding spin system which can be achieved via quantum annealing relying on the adiabatic theorem. Some of the inherent disadvantages of this procedure can be alleviated or resolved using a stochastic approach, and by coupling to the external environment. We show that careful engineering of the system-bath coupling at an individual spin level can further improve annealing.


Introduction
Quantum computation is expected to outstrip its classical counterpart in certain mathematical algorithms such as integer factorisation [1], in data searches [2], and also in large scale optimisation [3], especially if the target function corresponds directly to the quantum-mechanical description of the system, e.g., ground state [4] or correlation structure [5].Hence, should a challenging combinatorial problem be mapped onto a quantum system, a quantum simulator may indirectly solve the original problem in reasonable time.
A prototypical and suitable physical system is an Ising spin system, which suits experimental realisation, and all details of the hypothetical problem (e.g., the travelling salesman problem [6]) are encoded in the interaction strengths between spins, J ij .By finding the ground state, for given J ij , one solves the original problem [7][8][9][10][11][12][13][14].Of course, finding the ground state is itself an NP-hard problem for a large Ising system [15], in general.However, a potentially feasible strategy is to let Nature find the solution in a quick experimental protocol.
One approach-exploiting Kato's theorem [16]-is to prepare the system in the ground state of a simple Hamiltonian, and tune this Hamiltonian to the one whose ground state we are actually seeking.Kato's adiabatic theorem states that if a self-adjoint operator has an isolated eigenvalue with a potentially degenerate eigenspace and this eigenvalue does not split under a smooth self-adjoint perturbation, then there is a unitary transformation which transforms the unperturbed eigenspace to the corresponding perturbed eigenspace.In other words this theorem guarantees that a ground state of an unperturbed Hamiltonian can be evolved into the ground state of another Hamiltonian as long as the change of Hamiltonians happens smoothly enough.
One may apply a strong magnetic field in the x-direction to align all spins initially, then turning off the field slowly, and eventually measuring the spins in the z-directions.Unfortunately, in this process the energy levels of the instantaneous Hamiltonian usually become very close, and hence the driving must slow down to remain adiabatic (viz.critical slowing down [17][18][19][20][21]).Alternatively, proceeding in a finite time elevates the probability of the system being excited.As real quantum systems always interact with their environment, one may make a virtue of necessity, and allow the spins to couple to the thermal or quantum fluctuations of this external bath [22][23][24][25].Repeating the experiment a number of times yields an ensemble of final states from which we can pick the lowest energy configuration since calculating the energy for a given configuration is simple.In this approach, the likelihood of 'not finding the true ground state' diminishes exponentially with the number of trials [26].
Engineering the interaction with the environment can improve quantum annealing [27][28][29][30][31][32][33][34], and such active design, perhaps individually to each qubit, forms the heart of this paper.Inspired by energy transport optimisation through site-dependent coupling to the environment in photosynthetic systems [35,36], or individually addressable qubits on silicon [37][38][39], we propose the use of site-specific dissipation.We investigate this approach by evolving an Ising system via the Bloch-Redfield master equation.Within its approximations quantum annealing can achieve high efficiency [40][41][42][43][44][45][46].We show that online learning and adjustment of the individual coupling of each qubit to the environment indeed increases the probability of success, and that the annealing time is also an important factor in the protocol.

Model
Let us consider a collection of N identical, interacting qubits arranged at the vertices of an undirected graph, G, as depicted in Fig. 1: J ij σ z i σ z j with ferromagnetic couplings.A small field, m 0 = 8 × 10 -2 , is pinning a single qubit, σ 0 , in order to resolve the degeneracy of the ground state.The system is controlled by timedependent external magnetic field where τ is the annealing time, and m = μB 0 is chosen to overwhelm all other terms at the start of the experiment, i.e., B 0 max(J ij ).Here μ is the magnetic moment of a single qubit, and it is considered to be unity in the simulation, while the maximal magnetic field is chosen to be B 0 = 4.We consider two cases: (a) uniform coupling (J ij = J), and (b) random couplings, where every J ij are drawn from a normal distribution N (2, 0.1).We distinguish three classes of graphs: random non-complete graphs, a linear graph and a complete graph.The latter two are unique by their adjacency structure and form the extremes of connected graphs: a linear graph has the least, while the fully connected graph the most edges while being connected.In experiments we envisage a system somewhere in between these cases.In order to investigate these 'more typical' non-complete graphs, we sampled connected graphs for which half of all possible edges were missing.An exhaustive simulation could explore random connected graphs with k edges, but uniformly sampling this set of graphs is not a trivial task [47].
The Ising system is brought into contact with an infinitely large bath in thermal equilibrium modelled as The summation runs over the oscillator modes, while b † λ and b λ are the bosonic creation and annihilation operators for mode λ.Finally, the qubit-environment interaction is where we opted for A i = κ i (σ + i + σ - i ) with site-specific coupling κ i , and B = λ (b † λ + b λ ) as operators of the quantum system and of the environment, respectively.For the sake of simplicity, although unphysical, we assume κ i being independent of energy, i.e., a qubit is coupled to all bath modes equally.
We mention an alternative approach for regulate the system on a site-dependent basis, in which each qubit is controlled by its own transverse field [48], rather than an individually controlled interaction strength.Within a mean-field-like model Susa et al. showed that with such a spatio-temporal inhomogeneous control-field the detrimental first-order phase transition, present in a uniform system, can be avoided, hence the performance of the quantum annealing protocol can be improved.
The Hamiltonian, H = H G + H ext + H bath + H int , governs the time evolution ρtotal = -i[H int , ρ total ] in the interaction picture.As we are interested in the dynamics of the qubits, we trace over the bath states and obtain In solving Eq. ( 1) we make the three standard assumptions and follow Refs.[49][50][51].
The qubit-bath coupling is weak (Born approximation).The bath is initially uncorrelated.
The density matrix can be factorised (i.e., no appreciable correlation between system and bath) which is kept throughout the evolution, hence ρ total (t) ≈ ρ(t) ⊗ ρ bath , where ρ bath ∝ exp(-βH bath ) is the canonical density matrix at inverse temperature β.Finally, the bath is assumed to be Markovian, which allows us to change the integration limit in Eq. ( 1) from t to ∞ leading to Here C bath = Tr bath (B(t)B(t )ρ bath ) is the bath time-correlation function.Using the energy eigenbasis of H sys = H G + H ext , e.g., a|A(t)|b = A ab e iω ab t , we transform Eq. ( 2) into a matrix equation where is the noise power-spectrum of the bath.The function J(ω) = ηωe -ω/ω c is the Ohmic spectral density function with a cut-off frequency ω c and a dimensionless parameter η [52], while n(ω) = (e β ω -1) -1 is the mean oscillator number in mode ω at inverse temperature β.In the following we use κ(ω) = η/4πωe -ω/ω c (1 + n(ω)), thus the system couples differently to different modes of the bath.The aim of annealing is finding the ground state of a known Hamiltonian, hence in the following we report that value of κ (or κ = mean({κ i })) which it had at the moment when the probability of being in the instantaneous ground state stopped decreasing after the initial decline (see later in Figs. 2 and 3).Usually, for cold enough baths the relaxation rate is considered to be Ohmic [53], and, alternatively, an ensemble of coherent two-level systems can also reproduce Ohmic excitation spectrum [54,55].
We emphasise: the bath correlation time, τ bath , must be short enough, such that In other words, the fast oscillating quantities average out, and the terms with frequency ω abω cd will not give any significant contribution to the system evolution by t such that |ω abω cd | -1 t τ relax .This condition also sets a limit on the time-scale of the simulation [51,56], as beyond τ relax the observables do not change any more.
As a crude approximation τ relax ∝ κ -2 [57,58] while τ bath ∼ = max ( 2π ω c , β).Hence as long as τ bath τ relax is fulfilled, Eq. ( 3) should be an adequate description.The suitable β also depends on the graph and we use β = 22 for the linear graph, 30 for the random noncomplete graph, and 40 for the complete graph.These β values lead to κ i 0.007.

Results
Fidelity is a central measure in quantum computation and an upper bound can be calculated as [59,60] where ρ 1 and ρ 2 are any two density operators.In the following ρ 1 is fixed by the instantaneous ground state, while ρ 2 = ρ(t) is the instantaneous density operator.We also follow the time evolution of the instantaneous energy of the system, ε(t) = Tr(ρH sys (t)), compared to Tr(ρH sys (t = ∞)).
As a benchmark, we prepare the 7-qubit chain in its ground state with constant J ij = 1, and anneal with different τ without coupling to the environment.The probability of being in the instantaneous ground state, P, together with the logarithm of the gap between the ground state and the first excited state (dashed line), log( ), are shown in the Fig. 2. The gap initially diminishes rapidly, reaches its minimal value around t/τ ≈ 0.5, and then levels off at a value on the order of unity.Not surprisingly, around the min( ) fidelity drops as different states start mixing.However, for slower annealing (τ = 2, 4) scattering into the excited states is less, and higher P values are maintained, as expected from earlier studies [30,[61][62][63][64][65], hence τ plays an important role [57].Although short annealing time reduces Figure 2 The time dependence of P is plotted for an isolated 7-qubit chain (J ij = 1) for different τ (solid lines).
The dashed line is the logarithm of the energy gap, , as a function of scaled time Figure 3 (Top) Probability of being in the ground state (solid lines) over time for isolated (κ = 0) and open system (κ i ≡ κ = 0.001) for 7 qubits placed on an incomplete graph, where not all of them interacting with each other.The instantaneous energy (dashed lines) of both systems approaching the target energy.(Bottom) Populations of all 2 7 energy eigenstates at the end of the experiment for both systems probability of being in the ground state, we do aim to anneal fast, and hope that coupling to the environment helps the system to release its surplus energy and regain some of its lost population in the ground state.
Let us compare the population of the ground state for open and closed systems with fixed τ = 1  2 and uniform coupling to the environment.Figure 3 shows ε(t) and P for closed (κ = 0) and open (κ i ≡ κ = 0.001 for all i) systems.The open system coupled uniformly and weakly to the environment ends up in the instantaneous ground state with higher probability, and simultaneously the instantaneous energy, ε(t) approaches the target energy.The bottom panel of Fig. 3, compares population of states of each system: the ground state population remains much higher for a weakly coupled system than for the isolated system.Concluding, weak coupling to the environment may help the system to evolve into its ground state with higher probability.
Figure 4 shows P fin for fixed non-optimal coupling and for κ opt i .It is apparent that P fin < P optimised for a non-optimal κ, and also the final energy of the optimised case is closer to the target energy (the true ground state energy) than for the non-optimised system.
Next we vary κ i to maximise the population of the ground state at the end of the simulation, P fin .First, we consider a uniform coupling, κ i = κ opt , while later we allow for sitedependence, κ i = κ opt i .The evolution of P fin and the average coupling strengths, κ, are depicted in Fig. 5.We pick {κ i } corresponding to the highest P fin .Each simulation consists of two parts: one solves the Bloch-Redfield master equation and measure the fidelity between two given states, and the other optimises the target function y = 1 -max(fidelity) 2 with the constraint 0 < κ.We rely on the python3 module scipy.optimizeimplementing the Nelder-Mead algorithm [66].A run starts with an initial guess for κ, measures F as in Eq. ( 5), and then y is evaluated.
One may ponder upon the stability of the ground state population above some κ.We have simulated the time evolution for both non-optimised and optimised system-bath couplings.Figure 6 shows the results for a non-complete graph with 5 nodes coupled shown as a function of the average coupling strength κ.The graph contains N = 6 qubits, and it is a non-complete graph as described in the text.The graph is fixed for all simulations.The maximal value of P fin is reached around κ = 0.025.The red and blue dashed lines depict fin for the optimised and non-optimised cases, respectively.The gray dashed line shows the target energy.Here fin is calculated from the state when the annealing was stopped.(Bottom) The histograms illustrate the populations of states for the optimised and non-optimised (fixed κ) cases.The more pronounced bunching around the ground state, as expected, is visible for the optimised coupling strengths Figure 5 The probability of being in the ground state at the end of the simulation for a non-complete graph for random J ij ∼ N (2, 0.1).The coupling to the bath is either kept constant or it is optimised.The histograms show the populations of each energy eigenstates at the end of the evolution.Other parameters: N = 6, β = 30, and ω c = 30 to a bath at inverse temperature β with 0.001 ≤ κ i ≤ 0.007.If the annealing is stopped at t ∼ 3.2, one attains the highest χ = F 2 (ρ gs , ρ(t)) which we interpret as the probability of being in the instantaneous ground state.If the process continues, the system ends up with a smaller population in the ground state than expected from Boltzmann's distribution.While the optimised system undergoes a qualitatively similar evolution, χ reaches a higher value even though there is a period (see 1.5 t 2.7) when other couplings temporarily achieve higher χ values.
We have also analysed whether annealing is better compared to suddenly turning m(t) off at the very beginning of the time evolution after preparing the initial state.We chose a non-complete graph with 5 nodes, fixed J ij values, and varied τ from 0 to 2. The exper-Figure 6 Time evolution of χ for a non-complete random graph for 5 qubits and fixed τ = 0.5.All five qubits are coupled to the environment with the same strength, hence all κ i are equal to κ.The gray dashed line shows the probability expected for the final ground state in thermal equilibrium.The coloured solid lines are for non-optimized open systems with 0.001 ≤ κ ≤ 0.007.The top solid black line is for the optimised system while the bottom solid black line corresponds to an isolated system Figure 7 Comparison of annealing and non-annealing processes by varying the annealing time, τ .The solid lines are for open systems with κ = 0.003 (top) and 0.02 (bottom).The dash-dotted lines correspond to optimized κ, while the dashed lines are for an isolated system iments are repeated for two coupling strengths κ = 0.003 and 0.02.The results are plotted in Fig. 7 by line-triplets using the same colour for corresponding systems.One may draw three conclusions.First, all open annealed systems achieves higher χ values than the 'no-annealing' protocol, although for strong environmental coupling the excess probability is small.Secondly, a slow protocol (τ = 2) may outperform asymptotically all quicker quenches and optimisation does not improve χ appreciably.Finally, for all τ > 0 the time evolution of χ starts flat ( d dt χ ≈ 0), contrary to the 'no-annealing' protocol which starts increasing immediately.Let us now focus on the coupling of medium strength, κ = 0.003 (top panel of Fig. 7).Comparing lines with identical colours, one may observe that their order changes as τ varies.For a quick protocol, τ = 1  2 , the closed system performs the worst, an open system can improve on χ around t ≈ 3, while optimising κ further increases χ .For an intermediate value, τ = 1, the closed system outperforms the open system for long enough experiments, but it is still worse than the optimised system.Driving the magnetic field even slower (τ = 2) the closed system and the optimised system performs identically, i.e., optimisation does not improve the outcome.However, in reality no isolated system can be prepared, hence the conclusion remains: optimisation improves χ and may worth the effort.
The bottom panel of Fig. 7 corroborates our remarks.It is worth making two further comments.First, 'no-annealing' protocol seems to reach its asymptotic value very quickly for stronger couplings, hence repeating it twice or thrice, instead of any smooth quenches, may be a better approach.Second, the maxima of χ for open systems with annealing protocols are significantly reduced and they perform only slightly better than the 'no-annealing' protocol.
For both optimised and non-optimised cases we calculate the cumulative probability, P c , of 'not finding the true ground state in n consecutive experiments' as where P i is the probability of ending up in the ground state in the i th experiment.In experiments, where the system-bath coupling is not optimised for, P i remains constant throughout, hence P c (n) = (1 -P 1 ) n .However, with optimisation P i changes in each iteration and one may achieve a faster decrease than in the non-optimised case.Figure 8 demonstrates the difference between these approaches: in the optimised case the probability of missing the ground state tends to zero faster than the non-optimised case.As we learn better κ i values after each iteration, the probability of finding the ground state is higher than in the Figure 8 The cumulative probability, P c , is shown for four runs, with random initial couplings.In two runs (red) these initial couplings are frozen and kept constant, while in the other runs κ is optimized (blue) in each step.Even for the 'bad' initial guesses the curves decrease faster.(Inset) Same data with logarithmic ordinate non-optimised case.We note that, since the submission of our manuscript, some recent studies reached similar conclusion to ours, for example in terms of greedy optimisation of interaction parameters [67] or via genetic optimization of the annealing schedule [68].

Conclusion
We have studied the annealing of an archetypal transverse Ising spin (qubit) system coupled to an infinite heat bath with which it exchanges energy.We focused upon the probability of finding the ground state by the end of an experiment in finite time, if repeated runs are permitted.While these repetitions may extend the overall run-time of the annealing procedure, the probability of missing the true ground state in all experiments diminishes exponentially.We have also analysed the role of the annealing time, τ , and that of the site-dependent qubit-bath coupling strengths.However, it is left for future research to investigate the effects of graph theoretical quantities, such as adjacency structure, centrality, in-betweenness, etc, on the annealing procedure.
Focusing on a single experiment, we identified parameter ranges in which the environment can assist and improve the performance of the annealing.We have shown that the quench parameter, τ , and the system-environment couplings, {κ i }, can be optimised to improve the annealing process and keep the state of the quantum system close to the instantaneous ground state throughout the entire experiment.
Finally, repeating the simulation multiple times, the system ground state can be identified quicker if the system-bath coupling strength is varied in a supervised way at each iteration, compared to if one repeats the simulation with the system-bath coupling strength maintained constant.

Figure 1
Figure 1 Schematics of 7 qubits (black dots) interacting with each other with strengths J ij (olive solid lines) and with a Markovian bath (red springs) with couplings κ i and ω a , ω b are the eigenfrequencies corresponding the eigenstates |a and |b , respectively.We discard fast oscillatory terms and keep only those terms in the summation for which |ω abω cd | -1 τ relax .For clarity, we introduced the Bloch-Redfield tensor R abcd = -1 2 jk (δ bd r jk ac + δ ac r jk db ) with r jk ac = n A j an A k nc S(ω cn ) -A j ac A k db S(ω ca ) and r jk db

Figure 4 (
Figure 4 (Top) The probability of being in the instantaneous ground state at the end of the simulation, P fin , is