A thermodynamic approach to optimization in complex quantum systems

Optimization problems arise in different fields as computer science, physics, computational biology or drug design. One possible approach to solve an optimization problem with quantum algorithms is to encode the problem into an Ising Hamiltonian

$$\begin{equation} H_p = -\frac 1 2 \sum_{i,j} J_{ij} \sigma_i^z\sigma_j^z -\sum_i h_i \sigma_i^z, \end{equation} \tag{ 1 }$$

where $\sigma_i^w$ is the Pauli $w = x,y,z$ operator of the ith qubit, while the coupling strength J_ij and external field h_i specify the optimization problem. Models like the one in equation (1) were originally introduced in condensed matter physics to study disordered magnetic materials [1], but are now used to describe, e.g. neural networks [2], protein folding [3], socioeconomic behavior such as drug and tobacco use [4] and trading models in stock markets [5]. For reviews and perspectives on quantum computing applied to combinatorial optimization, see [6–8].

When the parameters J_ij and h_i in equation (1) are normally distributed, the model in (1) becomes a Sherrington–Kirkpatrick (SK) spin-glass [1, 9], whose energy minimization is an NP–hard problem. NP–hard problems are the hardest class of classical optimization problems [10, 11], examples include variants of the Traveling Salesman problem [12]. One disadvantage of NP–hardness as a concept is that it is a worst case statement, proven by mapping other NP–hard problems using polynomial resources. As a result, NP–hardness actually tells us nothing about the difficulty of a 'typical', in other words randomly generated, instance of a problem, and examples are known of NP–hard problems where typical instances are easy [13, 14]. For numerical studies such as this one, we also want the additional property of uniform hardness, which is a statement about the hardness of typical instances. The presence of a finite-temperature spin-glass transition [9] strongly suggests that as well as being NP-hard, the SK spin glass is uniformly hard, for this reason it has been used in other similar numerical studies [15, 16].

The use of quantum algorithms to solve such problems has the potential to provide a speedup over classical algorithms, provided that the system is made genuinely quantum by including a non-commuting driver term to the system Hamiltonian [17, 18]. The difficulty in proving speedups for NP–hard problems derives from the fact that the best classical algorithms are not known (strictly speaking it is not even known if the scaling is exponential for the best algorithm, although it is strongly suspected to be). In the case of unstructured search a more artificial problem where the best algorithm is known, Hamiltonians based adiabatic [19] and quantum walk [20] approaches (and in fact a family of algorithms which interpolate between the two [21]) are known to yield an optimal speedup equivalent to the one given by Grover's well known gate model algorithm [22, 23].

In order to find the ground state of the model described by equation (1), in this paper we use an approach that has been termed computing by cooling in [24]: in practice one employs a non-Markovian quantum bath in the form of a spin system with an easy-to-prepare low-temperature state, see a sketch of our proposal in figure 1. Coupling the system represented by the problem Hamiltonian (1) to such a cold bath, the system of interest is thus cooled down toward its ground state. Other thermodynamic-inspired approaches to classical and quantum computation have been previously proposed [25, 26].

Zoom In Zoom Out Reset image size

As a novel approach in the present paper we consider baths that can undergo quantum phase transitions (QPT) in the thermodynamic limit. We show that by tuning the parameter that determines the bath phase, we can speed up the search for the ground state of the problem Hamiltonian. Second–order QPT and other collective phenomena have indeed been shown to boost dynamic and thermodynamic performance in a number of thermal devices, both in classical and quantum regime [27–39].

While open system effects such as those experienced in our cooling approach could ruin the potential for quantum advantage, it is nevertheless known that even a projective measurement, which leads to pure decoherence, can lead to an optimal quantum speedup for unstructured search [40] (the problem encountered by Grover's algorithm) if performed in the correct basis. While unstructured search has the advantage of being a setting where quantum computing has known, provable speedup over classical computing, it has the disadvantage of behaving very different from the real optimization problem considered in [15], providing less information about possible performance on moderately sized real problems. For this reason, in this work, we elect to study a uniformly hard problem such as finding the minimum energy of SK spin glasses. The reason uniform hardness is important is that it implies that random instances will be hard. This is not true for typical NP-hardness, since the statement of NP-hardness is based on the ability to map to other problems and is thus a worst-case rather than a typical-case statement. There are in fact well known instances of NP-hard problems where randomly generated instances are easy to solve [13, 14]. Ising spin glasses however have a finite temperature spin glass transition, which strongly suggests that even random instances will be hard [9] (this property is known as uniform hardness). This subtlety had led to problems in previous works benchmarking quantum annealers, where it was realized that although solving Ising models when limited to the hardware graph of flux qubit devices is NP-hard, it is in fact not uniformly hard, and such problems are easy for classical Monte-Carlo algorithms [41].

This work is organized as follows. In section 2 we start by describing our model while in section 3 we discuss in detail the computational cooling scheme we propose. In section 4 we introduce a scheme which involves measuring the auxiliary system after its collision with the system of interest and in section 5 we summarize and discuss our findings.

In this paper we use the model described by the Hamiltonian in equation (1) with random coefficients as a test-bed to assess the performance of our algorithms and compare them to other quantum algorithms. In particular, we study a modified SK spin glass, as referenced in the introduction, consisting of randomly coupling all spins with coupling strengths and fields in equation (1) selected from a Gaussian distribution [42]. We start our analysis by following the procedure described in [43] as a basis for comparison. The parameters J_ij and h_i are drawn from a normal distribution with zero mean and variance equal to one. In particular, we complement the problem Hamiltonian (1) with a 'driver' Hamiltonian, and following previous works [15, 17, 18, 43], we will take it of the form

$$\begin{equation} H_d = -\sum_{i = 1}^{N} \sigma_i^x,\\ \end{equation} \tag{ 2 }$$

such that the total Hamiltonian reads

$$\begin{equation} H\left(\gamma\left(t\right)\right) = \gamma\left(t\right) H_d + H_p. \end{equation} \tag{ 3 }$$

Such a Hamiltonian results in quantum tunneling among the localized classical states, which correspond to the eigenstates $|E^{(i)}_p \rangle$ of H_p (forming the computational basis).

The ground state of H_d is the equal superposition of all the states forming the computational basis, a state that is in principle easy to prepare. Indeed, at t = 0, we prepare the system in the ground state of the driver Hamiltonian (2) and consider the following unitary dynamics generated by the Hamiltonian (3), a process sometimes called a 'two-stage quantum walk'. As time protocol we choose $\gamma(t) = \gamma_1 +(\gamma_2-\gamma_1) \theta(t-t_q)$ i.e. a quench in the driver's strength γ at time t_q , where we have used the step function $\theta(t)$, which is 0 for t < 0 and 1 for $t \unicode{x2A7E} 0$. The post-quench unitary dynamics of the total system with $H(\gamma_2)$ corresponds thus to a quantum random walk driven by H_d .

Notice that the system, initially in the ground state of H_d , effectively experiences a quench also at t = 0 when its Hamiltonian changes from H_d to $H(\gamma_1)$. We emphasize the importance of the two quenches, one at t = 0 and one at $t = t_q$, as means to change the system energy, which is otherwise conserved during the quantum walk dynamics.

We finally evaluate the expectation values of the different energy terms together with the fidelity with respect to the state of interest, as given by

$$\begin{equation} P\left(t\right) = \mathrm{tr} \left\{{\rho\left(t\right) \Pi^{\left(0\right)}_p}\right\}, \end{equation} \tag{ 4 }$$

with $\Pi^{(0)}_p = |E^{(0)}_{p} \rangle \langle E^{(0)}_{p} |$ the projector onto the problem Hamiltonian ground state.

The results for one specific realization of the disorder in (1) are shown in figure 2 where we plot the energy terms and the fidelity (4) as functions of the time. As a baseline comparison for the efficiency of the quench+walk process, we compare the energy and the fidelity to the corresponding values obtained for the ground state of H_p and of the total Hamiltonian $H(\gamma_2)$.

Zoom In Zoom Out Reset image size

In figure 3 we also show the results of an annealing process with a linear control function $\gamma(t) = \gamma_1+ (\gamma_2-\gamma_1) t/t_f$ for the same values of the control parameters γ₁ and γ₂ as in figure 3. Both figures 2 and 3 will be used as a basis of comparison for the cooling methods introduced and discussed in the following sections.

Zoom In Zoom Out Reset image size

We do not consider many of the factors which realistically exist in a large many-qubit annealing system. The reason we have elected to do this is so that our study can isolate the cooling effects from other competing effects which would greatly complicate our analysis. The largest consideration here is that the system will likely be in contact with an uncontrolled bath. These effects are important in real superconducting flux qubit architectures, the only fully controllable annealing systems built with thousands of qubits to date. These effects are not always detrimental, and benefits from interaction with a thermal bath have been experimentally demonstrated [44]. In fact, reverse annealing as implemented on D-Wave systems would not be able to solve optimization problems without dissipation [45, 46]. Realistically, bath effects will often still be undesirable; it has been demonstrated that coherent regimes can be reached experimentally by quenching very fast [47], but such approaches are unlikely to be practical in our setting since a large evolution time may be needed. Alternatively, in the future the protocols described here could be simulated on a fault tolerant universal quantum computer, thus avoiding having bath effects altogether. Simulations of annealing protocols on gate-model machines have recently gained interest through the approximate quantum annealing (AQA) algorithm proposed by Willsch et al [48].

In this section we introduce and discuss a computational cooling protocol to efficiently approach the ground state of H_p . To this purpose we let the system represented by the Hamiltonian H_p interact with an auxiliary system which is initially prepared in a state of low energy, see figure 1. Specifically we connect the spins of the SK model to an auxiliary 1D transverse-field spin-1/2 Ising chain, whose eigenvalues and eigenstates are well known [49]. The total Hamiltonian thus reads

$$\begin{align} H_{\mathrm{tot}}& = H_p+ \alpha H_A+H_I, \end{align} \tag{ 5 }$$

$$\begin{align} H_A& = -\sum_{i = 1}^{N} \left(1-f\,\right) \Sigma^x_{i}\Sigma^x_{i+1}-f \Sigma^z_{i}, \end{align} \tag{ 6 }$$

$$\begin{align} H_I& = -J\sum_{i = 1}^{N} \sigma^x_{i}\Sigma^x_{i}, \end{align} \tag{ 7 }$$

where we have introduced the Pauli operators $\Sigma^w_{i}$, $w = x,y,z$ for the auxiliary system, which we assume to be characterized by open boundary conditions ($N+1\equiv 1$ in equation (6)). The dimensionless parameter α is introduced to lower the value of the ground state energy of the auxiliary system, and as we will see, to enhance the cooling effect when α > 1.

It is worth noting that in the thermodynamic limit the auxiliary system in equation (6) exhibits a second order QPT at the critical point $f_c = 1/2$ (see [49]), with a doubly degenerate ground state and non-vanishing magnetization along the x-axis for $f\lt1/2$.

We initially prepare the system in the ground state of (2), and the auxiliary system in its own ground state $|E^{(0)}_A \rangle$. The two systems are then allowed to evolve following the unitary dynamics generated by the total Hamiltonian $H_\textrm{tot}$ for a time interval $\Delta t$, after which we disconnect the auxiliary system, and connect the system to a new auxiliary system freshly prepared in its ground state $|E^{(0)}_A \rangle$. We then iterate this procedure for a number n_c of cooling cycles. The resulting problem's energies of this cooling procedure, for different values of the parameter f can be seen in figure 4. First of all, we notice that for all f > 0 we observe some cooling effect in the system. A closer inspection of the results clearly indicates that the cooling procedure is most efficient for $f\simeq 0.6$, both in terms of cooling rate of $\langle {H_p} \rangle$ and of increasing value of the fidelity. While the critical value of f for the Hamiltonian H_A is $f_c = 1/2$, such a value is renormalized by the interaction with the system through the interaction Hamiltonian, and the value of f required to induce disorder in the auxiliary system becomes larger, $f\simeq0.6$, see the discussion in appendix A. We have checked that adding a term of the type $\sigma^y_i\Sigma^y_i$ to the Hamiltonian (7) does not lead to any substantial change to the results shown in figure 4.

Zoom In Zoom Out Reset image size

Comparing the results of figure 4 with those of figure 2, we see that that connecting the system to the external cooler improves significantly the optimization of the problem's energy. The lowest energy obtained with the repeated collisions are similar to that obtained with the annealing reported in figure 3 in a similar timescale.

We also evaluate the long time behavior of $\langle {H_p} \rangle$ as a function of f, as given by the final values of $\langle {H_p} \rangle$ at the end of the last cooling stroke. The results are shown in figure 5: inspection of this figure confirms that the maximal cooling effect is obtained for $f\simeq 0.6$.

Zoom In Zoom Out Reset image size

In appendix B, we complement our study by investigating the effect of the parameter α renormalizing the auxiliary system Hamiltonian (see equation (5)) on the cooling protocol. Furthermore, we consider the effect of freezing the system dynamics by reducing the interaction strength J appearing in equation (7). We also show that the long time behavior of the system undergoing repeated collision with the auxiliary bath is practically independent of the chosen initial state. Finally, in the same appendix, we consider different realizations of the quenched disorder in the SK chain (different sets of values of J_ij and h_i ): the cooling performance of our protocol does not show any substantial change, see figure 14.

In order to evaluate the purity of the system across different strokes, we compute the von Neumann entropy of the reduced density matrix of the system:

$$\begin{equation} S\left(\rho_S\right) = -{\mathrm{tr}} \rho_S \text{ln} \rho_S. \end{equation} \tag{ 8 }$$

In figure 6, we report $S(\rho_S)$ as a function of the time. After the first cycle the system gets close to the maximally mixed state, corresponding to a large system-environment entanglement. Subsequently, $S(\rho_S)$ decreases as the number of cooling cycles increases. Thus, while the repeated collision protocol reduces the energy of the system, such a cooling is achieved thanks to the establishment of quantum correlations with the auxiliary system. This was also noted in smaller collision models [50]. In the inset of the same figure we plot the fidelity as a function of the time for ten cooling cycles.

Zoom In Zoom Out Reset image size

It is worth noticing that the process describing the evolution of the system state is non–Markovian. Indeed we use a finite interaction time $\Delta t = 5$ in dimensionless units which does not lead to a Gorini–Kossakowski–Susarshan–Lindblad (GKSL) master equation for the system of interest. The evolution of an open quantum system is Markovian if and only if the process is described by a GKSL master equation, see, e.g. [51] and references therein. The GKSL master equation is retrieved if one considers vanishing interaction time $\Delta t$ between the bath and the system, and interaction strength that scales as $H_I\sim1/\sqrt{\Delta t}$ [52]. In appendix B we investigate the performance of the cooling protocol in that regime and find that the finite time interaction protocol performs considerably better in a single interaction stroke than the Markovian one with several short collisions. Thus we conclude that the bath and the system need to interact for a finite time for the bath to extract an appreciable amount of heat from the system.

Quantum thermalization machines made of ancillary spins were also proposed in [53, 54] to realize the thermal states of single spins or finite many-body systems. In those references, the authors also used repeated interactions with engineered environments, however a fine tuning of the baths energy level was required to achieve the target thermal states.

It is worth briefly discussing the practical implementation of the protocol depicted in figure 1. The most naive approach is to consider a new spin chain being coupled for each cooling cycle. Such an implementation would be impractical, since the number of qubits would grow linearly with the number of cooling cycles. Alternatively, if we had a separate procedure to 'reset' the same spin chain to an approximate ground state, we could reuse the same chain, see for instance [55]. This would introduce a period of time during which the chain was reset and the spin glass was allowed to undergo dynamics. Based on numerical observations in [15], it is likely that during this time period the system would rapidly equilibrate toward a long-time average, so adding such a time period would be unlikely to be detrimental. Adding such a time period would however complicate our model and give us another parameter to consider so it is therefore undesirable at the level of detail which we are working. We argue that the cooling time could be avoided using a constant number of spin chains using the following procedure: after each spin chain has interacted with the system it immediately begins the resetting procedure, meanwhile another already reset spin chain interacts with the system. As long as the number of spin chains is greater than the ratio of the resetting time to the time of each cycle, a reset spin chain will always be available, see also discussion in [56].

One question is how the coupling could be turned on or off in practice, but this will depend on the underlying system. For example, in circuit based flux qubit architectures, the coupling elements are the same as the qubits, but biased into a monostable rather than bistable regime [57]. In such an architecture it would be relatively straightforward to tune coupling by independently adjusting the biases on the different qubits (and single-body correction factors on the qubits they couple). In this setting a model of individually tuning couplers is realistic, although instantaneously turning off coupling may be difficult, making a model where the coupling is reduced non-instantaneously over time slightly more realistic. Simulating such systems would however be much more numerically demanding. Since we are interested in qualitative behavior of realistic systems rather than simulating real devices, instantaneous switching is sufficient for these purposes. Likewise, in atomic systems, coupling could be controlled by manipulating the physical separation between the system and bath qubits could be an effective method of turning coupling off, especially since all couplings are turned on and off at the same time. Shuttling of ions is considered to be an important component to trapped ion quantum computing and protocols to effectively and quickly move ions are being developed [58], so this could be a realistically achievable way of controlling the coupling.

As for the resetting procedure itself, we first note that, strictly speaking, the third law of thermodynamics forbids a system from being prepared exactly in its ground state in a finite time. For our purposes however, we assume that the spin chain is prepared in an approximate ground state with a finite energy gap, where, at a small but finite temperature, probabilities to be excited can be neglected. In fact a common approximation in quantum annealing is to assume that the system originates in the ground state of the driver Hamiltonian [59], in many ways our approximation is no different. As this state is an approximate thermal state, it is a free thermodynamic resource, see for instance [60]. From such a simple product state, it is possible to prepare strongly correlated systems that are ground states of quantum spin models, [61, 62]. This would require a number of gates that scales polynomially with the number of qubits in the environmental chain. This gates may require external thermodynamic work to be realized which therefore scales polynomially with the number of qubits.

We emphasize that our proposal can also be implemented in quantum circuit models as a quantum collision model, also known as a repeated interaction model. In fact, small quantum collision models have already been demonstrated in superconducting quantum circuits [55, 63]

In this section we consider again the system described in equations (5)–(7) and combine the collisional cooling method discussed in the previous section, with projective measurements on the auxiliary system. Let $|E_A^{(j)} \rangle$ be the eigenstates of H_A with projector $\Pi^{(j)}_{A} = |E_A^{(j)} \rangle\langle E_A^{(j)} |$, and ρ_tot the state of the combined system plus environment at a given time. The probability of measuring the energy $E_A^{(j)}$ is then given by

$$\begin{equation} p\left(E_A^{\left(j\right)}\right) = \mathrm{tr}\left[ \Pi^{\left(j\right)}_{A} \rho_{\mathrm{tot}}\right], \end{equation} \tag{ 9 }$$

and the post-measurement reduced state of the system is

$$\begin{equation} \rho_S\left(t|E_A^{\left(j\right)}\right) = \mathrm{tr}_A \left[ \Pi^{\left(j\right)}_{A}\rho_{\mathrm{tot}}\right]/p\left(E_A^{\left(j\right)}\right). \end{equation} \tag{ 10 }$$

One can then evaluate the expectation value of $\langle {H_p} \rangle$ conditioned by the outcome of measurement $E_A^{(j)}$ on the auxiliary system

$$\begin{equation} \langle {H_p} \rangle_{E_A^{\left(j\right)}} = \mathrm{tr}_S \left[\rho_S\left(t|E_A^{\left(j\right)}\right) H_p\right]. \end{equation} \tag{ 11 }$$

To evaluate the effect of a measurement after the repeated cooling strokes, we first apply the measure protocol to the cooling dynamics as depicted in figure 4: at the end of a given number of cooling cycles we evaluate the probability $p(E_A^{(j)})$ of measuring a given eigenstate of H_A , and the corresponding value of the post-measurement system energy $\langle {H_p} \rangle_{E_A^{(j)}}$. More precisely, we set $n_c = 1,\,2, \, 3$ or 4, cool the system as described in the previous section, and then measure the energy of A. The results of such a measurement protocol are shown in figure 7, for different values of the parameter f. We see that, in general, measuring H_A has no practical beneficial effect in the cooling of the system: the total probability that a measurement returns a value of $\langle {H_p} \rangle$ lower that the value at the end of the cooling stroke is smaller than 50%. However, at the end of the 4th stroke and for $\alpha = 2,\, 3$ there is a significant probability (${\gtrsim}{40}$%) that the measurement returns a value of $\langle {H_p} \rangle$ lower than the one at the end of the cooling stroke, resulting in additional net cooling.

Zoom In Zoom Out Reset image size

Next, we investigate whether the measurement of the auxiliary system can be used to speed-up the cooling of the system of interest and enhance the fidelity with respect to the target state equation (4) when compared to the simple collisional cooling shown in figure 4.

To this end, we implement the following cooling algorithm: at t = 0 we prepare both the system of interest and the auxiliary system as described above, i.e. the system starts in the ground state of (2) while the auxiliary cooler A starts in its own ground state. We then let the composite system evolve with the dynamics generated by the total Hamiltonian (5) for a time $\Delta t$, at the end of which we perform a measure of H_A , calculate the corresponding system post-measurement reduced state (10) and start a new cycle by connecting the system with a new auxiliary system freshly prepared in its ground state $|E^{(0)}_A \rangle$. We then iterate the above procedure for a given number of strokes. We consider two different protocols where we combine collisional cooling and measurements.

First, we only select post-measurement states corresponding to the second smallest eigenvalue of H_A , $E_A^{(1)}$, see figure 8. The rationale behind this choice is that, being the A system prepared in its ground state, finding it in its first excited state is the minimal requirement for the system of interest to be cooled. While this specific protocol might be of limited practical use, as the probability of achieving a sequence of exactly n_c measurements with $E_A = E_A^{(1)}$ can be pretty small, it is however a useful exemplification of the cooling mechanism obtained by energy transfer to the A system prepared in a cold state. This scenario might however be interesting in presence of a Maxwell daemon, selecting only the post-measurement state $| E_A^{(1)} \rangle$. The feasibility of this setup and its thermodynamic cost is however beyond the scope of the present paper.

Zoom In Zoom Out Reset image size

Second, we consider a more viable protocol consisting of a stochastic post-measurement selection protocol, more similar to a realistic experimental setup. At the end of each cooling stroke, we perform a measurement of H_A and as a consequence the system collapses in a post-measurement state $\rho_s(t|E_A^{(j)})$ with probability $p(E_A^{(j)})$, see equations (9) and (10). After such a projective measurement, we disconnect the system from the A system, and connect it to a new A system prepared in its ground state and iterate the procedure. The corresponding results are shown in figure 9. The curves in such a figure should be compared with those in figure 4 and 11 obtained with repeated collisions alone for the same value of the parameter f. Combining collisions and repeated measurements gives a similar decreasing rate for $\langle {H_p} \rangle$. Averaging over many stochastic 'trajectories' as the one depicted in figure 9, one would obtain the curves for $\langle {H_p} \rangle$ and P(t) shown in figure 4, which thus represent the expectation values of the quantities under scrutiny when repeating the stochastic post-measurement selection protocol several times.

Zoom In Zoom Out Reset image size

In this paper we have presented a thermodynamics-inspired protocol to find the minimum energy of quantum systems characterized by complex Hamiltonians. Specifically we have considered a classical optimization problem, that is turned into a quantum cooling process by using a non-Markovian bath that induces coherent cooling dynamics in the system of interest. The bath is initially prepared in its ground state and lowers the system energy through repeated collisions. The optimal working regime for our machinery is found in the range of parameters where the bath exhibits a QPT in the thermodynamic limit, thus emphasizing the importance of collective effects in thermal devices.

We conjecture that this result is not restricted to the specific bath considered in the present paper: any bath operating in its critical region should be extremely 'susceptible' to external disturbances and thus more capable of exchanging energy from the systems they are put in contact with.

By combining the cooling protocol with measurements on the bath alone, we are able to improve further the protocol, and in particular to keep the post-measurement state of the system pure. The possibility of using feedback control, mimicking a Maxwell daemon, with post-measurement selection rules affecting the bath alone seems to be an interesting direction for future explorations.

It would be interesting to compare to other algorithms, for example gate model variational algorithms like the quantum approximate optimization algorithm [16, 46]. Due to the similar continuous-time setting we have limited the present study to comparison with annealing and multi-stage quantum walks, but a broader comparison would be enlightening. A systematic analysis of the thermodynamic and computational resources required for our protocol and the comparison with alternative schemes remains an open problem.

While we studied the performance of the protocol in terms of cooling and fidelity rate for a wide range of parameters, it can be further improved by systematic optimization of the relevant parameters. In particular it might be interesting to consider baths with a number of spins different from the system of interest.

We acknowledge the support by the UK EPSRC EP/S02994X/1 and EP/T026715/1, and the Royal Society IEC\R2\ 222003.

A I gratefully acknowledges the financial support of The Faculty of Science and Technology at Aarhus University through a Sabbatical scholarship and the hospitality of the Quantum Technology group, at the School of Mathematics and Physics, Queen's University Belfast, during the initial stages of this project.

The data that support the findings of this study are openly available at the following URL/DOI: https://pure.qub.ac.uk/en/datasets/dataset-for-a-thermodynamic-approach-to-optimization-in-complex-q.

In the main text we have seen that the optimal cooling effect obtained with repeated collisions with the A system is obtained for $f\simeq0.6$, see equations (5)–(7). In this appendix we compare the behavior of the auxiliary system with (J ≠ 0 in equation (7)) or without interaction (J = 0) with the problem Hamiltonian. We remind the reader that the spectrum and the eigenstates of Hamiltonian (6) have been fully characterized in [49]. In particular, the ground state of H_A becomes doubly degenerate for f < 0 with non–vanishing magnetization along the x-axis for $f\lt f_c = 1/2$. In figure 10 we compare the expectation value of $\Sigma^x$ in presence or absence of interaction with the problem Hamiltonian for a system with N spins, together with the exact solution of the magnetization in the thermodynamic limit. We see clearly that the interaction with the problem Hamiltonian shifts the value of f for which the auxiliary system exhibits disorder toward a larger value than f_c . This must be compared again with the dynamical results of the collisional cooling, as exemplified in, e.g. figure 4 or 5, where we show very clearly that the optimal cooling effect is obtained for values of f larger than $f_c = 1/2$. Thus a comparison of figure 10 with figures 4 and 5 suggests that the optimal cooling regime occurs for values of f for which an order-disorder transition occurs in the auxiliary system coupled to the system of interest.

Zoom In Zoom Out Reset image size

In figure 11 we study the effect of the parameter α renormalizing the auxiliary system Hamiltonian (see equation (5)) on the cooling protocol. We see that increasing α improves the performance of the protocol introduced in this section, both in terms of the lowest problem's energy reached and the time to reach it. This is due to the fact that, for larger α, the system sees the auxiliary bath as colder, given that we initially prepare it in its ground state $|E^{(0)}_A \rangle$ as discussed above.

Zoom In Zoom Out Reset image size

In figure 12 we study the effect of freezing the system dynamics by reducing the interaction with the auxiliary system: specifically we implement two subsequent quenches on J, i.e. the interaction strength appearing in equation (7). Inspection of figure 12 confirms that decreasing the interaction strength between the system and the auxiliary cooler simply freezes the system in its state.

Zoom In Zoom Out Reset image size

In order to check that the repeated collision process works regardless of the initial state, in figure 13 we study the case in which the system to be cooled is initially prepared in its own ground state. From the figure, we can see that, after a few cycles, the results for this case and for the case analyzed in figure 11 approach each other and reach approximately the same behavior within the timescales depicted in the figure. Only for the large value f = 0.7, we observe some deviation possibly due to finite sizes. This strongly suggests that the system is reaching some kind of thermal equilibrium, independent of its starting state. This behavior is fundamentally different from unitary algorithms such as the adiabatic algorithm, where by construction memory of the starting state is maintained throughout the algorithm. It also suggests that the values here are the ultimate limit which the algorithm can achieve and that waiting for longer would not lead to improved performance.

Zoom In Zoom Out Reset image size

In figure 14 we show the result for the cooling protocol for three different realization of the SK chain, i.e. for different sets of values of J_ij and h_i which are drawn from a normal distribution with zero mean and unitary variance. One of the three sets is the one used in the main text.

Zoom In Zoom Out Reset image size

In figure 15 we show the difference in cooling performance of the protocol with long ($\Delta t = 5$) and short ($\Delta t = 0.1$) bath-system interaction strokes. This elucidates the need of non-Markovian memory effects in the environment. For the long-stroke protocol the prefactor of the interaction Hamiltonian (7) takes the value α = 3 as in the main text. In this regime, the environment interacts for a sufficiently long time such that back-flow of information from the environment to the system, a hallmark of non Markovianity, is possible. For the short-stroke protocol, the prefactor of the interaction Hamiltonian reads $\alpha = 3/\sqrt{\Delta t}$. In this regime the protocol dynamics for the system becomes equivalent to the one described by a Markovian GKSL master equation [52]. We see that the long-stroke protocol performs considerably better than the short one, indicating that bath and the system need to interact for a finite time for the bath to extract an appreciable amount of heat from the system.

Zoom In Zoom Out Reset image size