A thermodynamic approach to optimization in complex quantum systems

We consider the problem of finding the energy minimum of a complex quantum Hamiltonian by employing a non-Markovian bath prepared in a low energy state. The energy minimization problem is thus turned into a thermodynamic cooling protocol in which we repeatedly put the system of interest in contact with a colder auxiliary system. By tuning the internal parameters of the bath, we show that the optimal cooling is obtained in a regime where the bath exhibits a quantum phase transition in the thermodynamic limit. This result highlights the importance of collective effects in thermodynamic devices. We furthermore introduce a two-step protocol that combines the interaction with the bath with a measure of its energy. While this protocol does not destroy coherence in the system of interest, we show that it can further enhance the cooling effect.


I. INTRODUCTION
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 where σ w i 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 Eq. (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 optimisation, see [6][7][8].
When the parameters J ij and h i in Eq. ( 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 optimisation problems [10,11], examples include variants of the Traveling Salesman problem [12].One disadvantage of NPhardness 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 Eq. (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 Fig. 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 thermodynamicinspired approaches to classical and quantum computation have been previously proposed [25,26].
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

H S H A H I
FIG. 1: Schematic of the collision protocol: a series of auxiliary spin chains (squares) interact with the SK spin glass system (circles).Black and red solid lines connecting spins in the glass represent positive and negative pairwise interactions, respectively, as described by the terms proportional to Jij in Eq. ( 1).Their thickness represents their random relative magnitude.The gray broken lines represent the interaction between the spins in the system and in the bath, as given by J in Eq. (7).At the end of every collisional stroke the exhaust auxiliary chain is wasted and replaced by a new chain.
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 optimisation problem considered in Ref. [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 [41] (this property is known as uniform hardness).This subtlety had led to problems in previous works benchmarking quantum annealers, where it was realised 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 [42].
This work is organised as follows.In Sec.II we start by describing our model while in Sec.III we discuss in detail the computational cooling scheme we propose.In Sec.IV we introduce a scheme which involves measuring the auxiliary system after its collision with the system of interest and in Sec.V we summarise and discuss our findings.

II. PRELIMINARIES
In this paper we use the model described by the Hamiltonian in Eq. ( 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 Sherrington-Kirkpatrick spin glass, as referenced in the introduction, consisting of randomly coupling all spins with coupling strengths and fields in Eq. ( 1) selected from a Gaussian distribution [64].We start our analysis by following the procedure described in Ref. [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 such that the total Hamiltonian reads Such a Hamiltonian results in quantum tunneling among the localized classical states, which correspond to the eigenstates p ⟩ 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 γ(t) = γ 1 + (γ 2 − γ 1 )θ(t − t q ) i.e. a quench in the driver's strength γ at time t q , where we have used the step function θ(t), which is 0 for t < 0 and 1 for t ≥ 0. The post-quench unitary dynamics of the total system with H(γ 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(γ 1 ).We emphasise 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 (0) p | the projector onto the problem Hamiltonian ground state.
The results for one specific realization of the disorder in (1) are shown in Fig. 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(γ 2 ).
In Fig. 3 we also show the results of an annealing process with a linear control function γ(t) = γ 1 +(γ 2 −γ 1 )t/t f for the same values of the control parameters γ 1 and γ 2 as in Fig. 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.
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 optimisation 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].

III. COMPUTATIONAL COOLING
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 H(γ 2 ) , of the total Hamiltonian H(γ2), Eq. (3).Bottom: fidelity (4) as a function of time.We have used N = 9 spins.system which is initially prepared in a state of low energy, see Fig. 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 where we have introduced the Pauli operators Σ w i , w = x, y, z for the auxiliary system, which we assume to be characterised by open boundary conditions (N + 1 ≡ 1 in Eq. ( 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 Eq. ( 6) exhibits a second order QPT at the critical point f c = 1/2 (see Ref. [49]), with a doubly degenerate ground state and non-vanishing magnetization along the x-axis for f < 1/2.
We initially prepare the system in the ground state of (2), and the auxiliary system in its own ground state |E (0) A ⟩.The two systems are then allowed to evolve following the unitary dynamics generated by the total Hamiltonian H tot for a time interval ∆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 ⟩. 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 Fig. 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 ≃ 0.6, both in terms of cooling rate of ⟨H p ⟩ 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 ≃ 0.6, see the discussion in Appendix A. We have checked that adding a term of the type σ y i Σ y i to the Hamiltonian (7) does not lead to any substantial change to the results shown in fig. 4.
Comparing the results of Fig. 4 with those of Fig. 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 Fig. 3 in a similar timescale.
We also evaluate the long time behaviour of ⟨H p ⟩ as a function of f , as given by the final values of ⟨H p ⟩ at the

FIG. 4: Expectation value of the problem Hamiltonian (top)
and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system (6) for different values of f , with α = 3, nc = 5 cooling strokes of duration ∆t = 5 in dimensionless time units.The system is composed of N = 9 spins, and the parameters of Hp are the same as in Fig. 2.
end of the last cooling stroke.The results are shown in Fig. 5: inspection of this figure confirms that the maximal cooling effect is obtained for f ≃ 0.6.
In Appendix B, we complement our study by investigating the effect of the parameter α renormalizing the auxiliary system Hamiltonian (see Eq. ( 5) on the cooling protocol.Furthermore, we consider the effect of freezing the system dynamics by reducing the interaction strength J appearing in Eq. (7).We also show that the long time behaviour 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 Fig. 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: In Fig. 6, we report S(ρ S ) as a function of the time.After the first cycle the system gets close to the maximally mixed state, corresponding to a large systemenvironment entanglement.Subsequently, S(ρ 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.It is worth noticing that the process describing the evolution of the system state is non-Markovian.Indeed we use a finite interaction time ∆t = 5 in dimensionless units which does not lead to a Gorini-Kossakowski-Susarshan-Lindblad (GKSL) master equation for the sys-tem 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 ∆t between the bath and the system, and interaction strength that scales as H I ∼ 1/ √ ∆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 realise 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 Fig. 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 Ref. [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 Ref. [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 singlebody 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 behaviour 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 Ref. [60].From such a simple product state, it is possible to prepare strongly correlated systems that are ground states of quantum spin models, Ref. [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 realised which therefore scales polynomially with the number of qubits.
We emphasise 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]

IV. COOLING BY MEASURING THE AUXILIARY SYSTEM
In this section we consider again the system described in Eqs. ( 5)-( 7) and combine the collisional cooling method discussed in the previous section, with projective measurements on the auxiliary system.Let |E (j) A ⟩ be the eigenstates of H A with projector Π (j) A |, and ρ tot the state of the combined system plus environment at a given time.The probability of measuring the energy E (j) A is then given by p(E (j) and the post-measurement reduced state of the system is ρ S (t|E One can then evaluate the expectation value of ⟨H p ⟩ conditioned by the outcome of measurement E (j) A on the auxiliary system 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 Fig. 4: at the end of a given number of cooling cycles we evaluate the probability p(E (j) A ) of measuring a given eigenstate of H A , and the corresponding value of the postmeasurement system energy ⟨H p ⟩ E (j)

A
. 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 Fig. 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 ⟨H p ⟩ 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 α = 2, 3 there is a significant probability (≳ 40%) that the measurement returns a value of ⟨H p ⟩ lower than the one at the end of the cooling stroke, resulting in additional net cooling.
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 Eq. ( 4) when compared to the simple collisional cooling shown in Fig. 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 ∆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 ⟩. 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 , see Fig. 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 A 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 ⟩.The feasibility of this setup and its thermodynamic cost is however beyond the scope of the present paper.
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 postmeasurement state ρ s (t|E (j) A ) with probability p(E (j) A ), see Eqs. ( 9)- (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 Fig. 9.The curves in such a figure should be compared with those in Fig. 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 ⟨H p ⟩. Averaging over many stochastic "trajectories" as the one depicted in Fig. 9, one would obtain the curves for ⟨H p ⟩ and P (t) shown in Fig. 4, which thus represent the expectation values of the quantities under scrutiny when repeating the stochastic post-measurement selection protocol several times.

V. CONCLUSIONS
In this paper we have presented a thermodynamicsinspired protocol to find the minimum energy of quantum  5)-( 7)).In this protocol we select the post-measurement state corresponding to the second smallest eigenvalue of HA: E A .Top panel: expectation value of the system Hamiltonian Hp after the measurement as a function of the number of cooling cycles.Inset: probability of the selected state of the auxiliary system pA(E 1 A ). Bottom panel: expectation value ⟨Hp⟩ as a function of t.Inset: fidelity as a function of t. 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 postmeasurement 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.
FIG. 10: x-magnetization of the A system in presence or absence of the interaction with the problem Hamiltonian.The full line is the exact solution for ⟨Σ x ⟩ in the thermodynamic limit [49].
colder, given that we initially prepare it in its ground state |E A ⟩ as discussed above.In Fig. 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 Eq. (7).Inspection of Fig. 12 confirms that decreasing the interaction strength between the system and the auxiliary cooler simply freezes the system in its state.
In order to check that the repeated collision process works regardless of the initial state, in Fig. 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 analysed in Fig. 11 approach each other and reach approximately the same behaviour 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 behaviour 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.
In Fig. 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.
In Fig. 15 we show the difference in cooling performance of the protocol with long (∆t = 5) and short (∆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

FIG. 11: Expectation value of the problem Hamiltonian (top)
and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system (6) for different values of f and α, with nc = 5 cooling strokes of duration ∆t = 5 in dimensionless time units.The system is composed of N = 9 spins, and the parameters of Hp are the same as in Fig. 2. Larger values of α, e.g.α = 4, results in a deterioration of the cooling performance (data not shown).
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 shortstroke protocol, the prefactor of the interaction Hamiltonian reads α = 3/ √ ∆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.Expectation value of the problem Hamiltonian (top) and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system (6) for different values of f , with α = 3, nc = 12 cooling strokes of duration ∆t = 5 in dimensionless time units.The system is composed of N = 9 spins, and the parameters of Hp are the same as in Fig. 2. The dashed lines correspond to a time protocol where we perform a quench on J every four strokes, the continuous line correspond to constant interaction strength J = 1.The time protocol J(t) is plotted in the inset of the first figure.and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system (6) for different values of f , with α = 3, nc = 10 cooling strokes of duration ∆t = 5 in dimensionless time units.The full symbols correspond to the case where the system is initially prepared in the ground state of (2), as in figs.4-12, while the empty symbols correspond to the case where the system is initially prepared in its own ground state.

4 P
FIG. 2: Top: Expectation values of the different contributions to the total Hamiltonian (3) as a function of the time, with the quench+walk protocol, and with the ergotropy extraction protocol.The values of the control function parameters are taken to be γ1 = 4 and γ2 = 1.For reference, we also plot the ground state energy E (0) p of Hp, Eq. (1) and that, E

FIG. 5 :
FIG.5: Final value of ⟨Hp⟩ as a function of f .For each value of f , ⟨Hp⟩ fin is the final value of the system energy at the end of the last cooling stroke of Fig.4.

FIG. 6 :
FIG.6: Main panel: von Neunmann entropy of the reduced density matrix of the system S(ρS) as a function of the time for three values of f .The dashed line corresponds to the entropy of the completely random state N ln 2. Inset: fidelity as a function of the time for the same values of f .α = 3, nc = 10 cooling strokes of duration ∆t = 5 in dimensionless time units.The system is composed of N = 9 spins, and the parameters of Hp are the same as in Fig.2.

FIG. 7 :
FIG.7: Post-measurement system energy(11) as a function of pA as obtained by measuring HA at the end of the first four cooling strokes, for different values of f , with α = 3 and N = 9.The rest of are as in Fig.4.Dashed lines: ⟨Hp⟩ at the end of the cycle, i.e. immediately before the measurement, see also Fig.4for the cycles with α = 3.The percentage gives the fraction of measurements that return a value of ⟨Hp⟩ below the value at the end of the cycle.

FIG. 8 :
FIG.8: Collisional cooling strokes followed by a projective measurement of the auxiliary system Hamiltonian HA with α = 1, ∆t = 5 and = 1 (see Eqs. (5)-(7)).In this protocol we select the post-measurement state corresponding to the second smallest eigenvalue of HA: E

FIG. 9 :
FIG. 9: Collisional cooling strokes followed by a projective measurement of the auxiliary system Hamiltonian HA (see eqs. (5)-(7)).In this protocol we randomly select the postmeasurement state with probability p(E (j) A ), Eq. (9).Top figure, main panel: expectation value ⟨Hp⟩ as a function of t for a single realization of the repeated measurement process.Inset: fidelity as a function of t.Bottom figure: ⟨Hp⟩ as a function of t for 2 different values of α and f = 0.6.

30 P
FIG. 12:Expectation value of the problem Hamiltonian (top) and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system(6) for different values of f , with α = 3, nc = 12 cooling strokes of duration ∆t = 5 in dimensionless time units.The system is composed of N = 9 spins, and the parameters of Hp are the same as in Fig.2.The dashed lines correspond to a time protocol where we perform a quench on J every four strokes, the continuous line correspond to constant interaction strength J = 1.The time protocol J(t) is plotted in the inset of the first figure.

7 FIG. 13 :
FIG. 13: Expectation value of the problem Hamiltonian (top)and fidelity (bottom) as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system(6) for different values of f , with α = 3, nc = 10 cooling strokes of duration ∆t = 5 in dimensionless time units.The full symbols correspond to the case where the system is initially prepared in the ground state of (2), as in figs.4-12, while the empty symbols correspond to the case where the system is initially prepared in its own ground state.

7 FIG. 14 : 1 FIG. 15 :
FIG. 14: Expectation value of the problem Hamiltonian divided by the absolute value of the ground state energy as a function of t for the cooling protocol consisting of repeated collisions with the auxiliary system (6) for different values of f , with α = 3, nc = 5 cooling strokes of duration ∆t = 5 in dimensionless time units.The different curves with the same color correspond to different realization of the SK model, i.e. to different sets of Jij and hi, sampled from the same normal distributions.The full symbols correspond to the SK model used in the main text, the empty symbols correspond to two different realizations of the quenched disorder.