Benchmark of quantum algorithms for ground state preparation in the presence of noise

1 Introduction

Preparing ground states of many-body Hamiltonians is a central task in quantum simulation and computation. Various strategies have been developed, including adiabatic evolution [3], variational algorithms, such as the Quantum Approximate Optimization Algorithm (QAOA) [2, 4, 5], and dissipative or cooling protocols [1, 6–21]. In current experimental settings, all these approaches are limited by the presence of noise, including decoherence and imperfections, which degrades the quality of the prepared states and imposes limits on the reachable energies. Thus, the following questions naturally arise: how does the achievable energy scale with the noise rate, and which strategy performs best for given hardware constraints?

The adiabatic algorithm [3] prepares the target ground state by slowly evolving the system from the ground state of a simple initial Hamiltonian to the target Hamiltonian, relying on the adiabatic theorem to keep the state close to the instantaneous ground state throughout the evolution. Its performance is often limited by the minimum spectral gap along the interpolation path, which closes at quantum phase transitions and demands impractically long runtimes [22, 23]. Variational algorithms such as the Quantum Approximate Optimization Algorithm (QAOA) [2, 5] circumvent some of those issues by optimizing the parameters of a fixed-depth quantum circuit in order to minimize the energy expectation value. Their scalability is often limited by barren plateaus at large circuit depth [5, 24], or by the expressivity of the chosen family of gates. Cooling algorithms take a qualitatively different approach: they engineer a coupling between the system and a cold auxiliary bath and extract energy through repeated cycles of joint evolution and bath reset [9]. Closely related are dissipative state preparation (DSP) protocols [6–8, 25], which engineer the system-bath coupling so that the dissipative dynamics drives the system to a desired target state. Both cooling and dissipative algorithms are described in terms of a channel which is repeatedly applied. Their performance is limited, among other factors, by the spectral gap of that channel.

The last two algorithms mentioned above have the advantage that the state of interest is prepared as the stationary state of the evolution (i.e. the fixed point of a linear map). Several complementary approaches have been put forward in that context. For instance, dissipative state preparation (DSP) schemes have been proposed to prepare ground states of frustration-free Hamiltonians [6, 7]. Ding et al. and Zhan et al. [12, 26] proposed a filtered Lindbladian approach to prepare the ground or the thermal state of arbitrary (local) Hamiltonians (see also [20]). Lloyd et al. [10, 11] introduced both a cooling scheme that uses time-dependent system-bath couplings and a quasiparticle cooling framework underlying a recent quantum processor experiment [14]. In [1] a cooling protocol which uses static couplings that are instantaneously switched on and off, combined with randomized cycle times and multiple bath frequencies was introduced.

In this work, we analyze how noise affects some of the ground state preparation algorithms. We compare various algorithms in the presence of noise, with the goal of identifying the fundamental mechanisms that limit their performance and determining their optimal behavior. Our approach combines analytical derivations, where we obtain the achievable energy as a function of the noise rate $\kappa$, with numerical optimization of the control parameters for each algorithm. Prior work has studied noise effects on individual algorithms, including the robustness of adiabatic computation to decoherence [27–29] and noise-induced errors in quantum simulation. However, those works do not directly compare the performance of different algorithms.

We compare the plain adiabatic algorithm, QAOA and the cooling protocol introduced in [1]. As in that work, we choose a simple Hamiltonian that contains two different phases, one trivial and the other one topological, connected by a quantum phase transition [30–32]. This allows us to derive analytical expressions and to perform numerical simulations in sufficiently large system sizes. The initial state of all the protocols is a trivial one (product state, or a more realistic state close to the trivial one, see below) and we choose as a figure of merit the relative energy, which, as opposed to the global fidelity, is an intensive quantity which is well suited in the context of quantum simulation.

We divide our analysis in two parts. The first one compares analytically the adiabatic algorithm with the cooling protocol in the weak coupling regime, where one can derive an effective master equation describing the dynamics. In that case, we optimize the evolution time in both methods to reach the lowest energy. As expected, the relative performance of the methods depends on the phase of the Hamiltonian. We show that both the energy scaling with the noise rate $\kappa$ and the prefactor in front of it are relevant for analyzing each method.

In the second part we compare two algorithms in which several control parameters are optimized: QAOA and the optimizable cooling algorithm of [1]. Note that, thanks to the variational parameters, both outperform the adiabatic algorithm and the multifrequency cooling algorithm explained previously. Thus, we only concentrate on the optimized QAOA and cooling algorithm. In order to enable a fairer comparison, we optimize in both the same number of parameters. We find that both algorithms perform comparably in the trivial phase but cooling outperforms QAOA in the topological regime. Additionally, we observe that the cooling algorithm is more resilient to parameter imperfections. We emphasize that this paper only compares the final performance of the protocols, since our aim is to understand the fundamental limits of the different algorithms.

The remainder of the paper is organized as follows. Section 2 provides an overview of the main results and their physical interpretation. Section 3 defines the model, the bath, and the noise channel. Section 4 describes the three algorithms we consider here. Section 5 derives the analytic scaling laws with noise strength for cooling and adiabatic evolution and compares them across the phase diagram. Section 6 presents a numerical comparison of QAOA and optimized cooling, including an analysis of robustness to parameter imperfections. Section 7 summarizes our conclusions and outlook.

2 Summary of results

In this section we present and discuss the main results derived in this work. We briefly introduce the model Hamiltonian, the noise model, and the figure of merit, as well as the quantum algorithms we are going to analyze. We then consider the weak coupling regime and compare the cooling and the adiabatic protocols. Then, we relax the condition of weak coupling and present the results for QAOA and the cooling methods.

2.1 Model

The goal is to prepare the ground state of a family of translationally invariant quadratic fermionic Hamiltonians, $H_S(𝜃)$, parametrized by $𝜃\; \in \; [0,\pi ∕2]$. This is the same textbook model as the one used in ref. [1]. Due to translational invariance, $H_S(𝜃)$ can be analytically block-diagonalized via a Fourier transformation, leading to a direct sum of Hamiltonians acting only on the mode pairs $(k,-k)$ [32]. Each of these conserves the fermionic parity, which allows us to use a tensor product formulation of the fermionic states. In particular, we will call ${|0⟩}^{\otimes N}$ the vacuum state. For $𝜃\; =\; \pi ∕2$, the Hamiltonian is local and the ground state is the vacuum. For $𝜃\; =\; 0$ the Hamiltonian is at the topological point of the Majorana chain with zero correlation length [33, 34]. The model possesses a phase transition at $𝜃\; =\; \pi ∕4$ that separates the trivial ($𝜃\; >\; \pi ∕4$) and the topological () phases. At the critical point, the gap for the mode with $k\; =\; N∕2$ closes.

As a figure of merit for the different algorithms, we consider the relative energy

$$e\; =\left|\frac{E\; -E_{GS}}{E_{GS}}\right|,\text{ with }E\; =\; tr(H_S\rho )$$

(2.1)

where $E_{GS}$ denotes the ground-state energy of $H_S$ and $\rho$ the output state of the corresponding algorithm.

As in [1], the noise is modeled through a depolarizing master equation given in section 3.3, characterized by a noise rate $\kappa$. Note that it commutes with the dynamics generated by the Hamiltonian, which simplifies the analysis.

As mentioned in the introduction, we will consider three quantum algorithms: adiabatic, QAOA and cooling. There are many different ways of performing each of those algorithms, and it is out of the scope of this work to analyze every one of them. This is why we take some of the standard approaches used in the adiabatic and QAOA algorithms. In the first one, we start with the Hamiltonian at $H(\pi ∕2)$ and transform it linearly with time to $H(𝜃)$. In the QAOA, we consider pairs of layers that contain the evolution dictated by $H(\pi ∕2)$ and $H(0)$, and optimize with respect to the evolution times in each layer independently. Regarding the cooling algorithm, we select the one considered in our previous work [1] with a slight modification. Despite the fact that we consider a specific cooling algorithm, we expect that other cooling and dissipative algorithms will share many of the features that we obtain in our analysis.

A key distinction between the algorithms is the role of the initial state. Whereas in a cooling protocol the state of interest is the stationary state, in the adiabatic and QAOA approaches the relevant state is obtained via a unitary evolution. Consequently, the initial state plays an important role in the latter two algorithms, while in a cooling process the final state is essentially independent of it. As for the initial state, one could simply consider the vacuum state. However, this is already the ground state for $𝜃\; =\; \pi ∕2$, so that in that case both the adiabatic and the QAOA algorithm will give the perfect state (without performing any operation), whereas the cooling algorithm will give an imperfect one since it will have to reach a steady state. To enable a meaningful comparison in the presence of noise, it is crucial to adopt a common framework in which all algorithms are evaluated under the same noise model and their performance is optimized with respect to their respective control parameters. In particular, this requires considering realistic initial states that can be prepared under this noise, rather than idealized states such as ${|0⟩}^{\otimes N}$, thereby ensuring that the comparison reflects experimentally accessible conditions. Due to that, we use as initial state ${\rho }^{cool}(\pi ∕2)$, obtained from a simple cooling protocol at $𝜃\; =\; \pi ∕2$. Although this does not change the scalings that we find, it may slightly influence the numerical results.

2.2 Comparison between adiabatic and cooling algorithms

To understand how the adiabatic and cooling algorithms compare in the presence of noise, it is useful to analyze both methods analytically. Our central analytical result is the derivation of scaling laws for the achievable relative energy as a function of the noise rate $\kappa$. We are interested in the limit of low noise rate, and in the energy scaling with $\kappa$ in that limit. When comparing two algorithms, we will say that one has better scaling with $\kappa$ if the relative energy increases more slowly in terms of that parameter. For instance, a scaling with ${\kappa }^{1∕3}$ is better than one with ${\kappa }^{2∕3}$, despite the fact that $\kappa \; \ll \; 1$. As we show, this scaling depends on both the algorithm under consideration (adiabatic versus cooling) and on the phase being studied (trivial versus topological).

Cooling. We consider the multi-frequency, randomized-time cooling protocol developed in ref. [1] (see also section 4.3 and appendix A), where an external bath mode is locally weakly coupled to the system with a coupling constant $g$ and extracts energy iteratively. Each cooling cycle consists of coupling the system to a bath mode of transition energy $\mathrm{\Delta }$ for an interaction duration $\tau$, followed by resetting the bath. We analyze the state after the concatenation of $\ell$ such cooling cycles, where in each cycle the frequency $\mathrm{\Delta }$ and interaction duration $\tau$ are chosen randomly from specified distributions. We consider a uniform distribution for $\tau \; \in \; [0,2t]$, whereas for $\mathrm{\Delta }$ we take both continuous and discrete distributions in the interval $\mathrm{\Delta }\in \; [{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$, which contains the whole single-particle spectrum of energies $\{{𝜖}_m,\dots ,{𝜖}_M\}$. We further constrain our analysis to the weak-coupling (${(gt)}^2\ll \; 1$), cooling (${𝜖}_{m}t\; \gg \; 1$), and low-noise ($\kappa t\; \ll \; 1$) limits. Note that this requires $\kappa \; \ll {𝜖}_m$. As shown in appendix A, in the limit of $\ell \; \to \; \infty$, the state obtained by such a cooling process is close to the stationary state of an averaged map for any distribution of $\tau$ and $\mathrm{\Delta }$.

Let us briefly comment on the choice of distributions of the bath frequencies considered here. As shown in [1], for a single bath frequency $\mathrm{\Delta }$ (i.e., the distribution of the bath frequency is a delta function), the cooling rate as a function of the energy is described by a Lorentzian function proportional to $t$, centered at $\mathrm{\Delta }$, and with a width proportional to $1∕t$. Thus, only the modes with energies close to that frequency $\mathrm{\Delta }$ will be effectively cooled. As $t$ increases, the energy of the corresponding mode becomes lower, but fewer and fewer modes are cooled. This is why it is necessary to choose many frequencies $\mathrm{\Delta }$ in an interval that includes all mode energies. One possibility is to consider a continuous distribution $\mathrm{\Delta }\in \; [{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$. This is analyzed in appendix A and allows us to obtain analytical expressions for the cooling rate and the final energy. Alternatively, one can consider a distribution containing $R$ discrete frequencies, equispaced in $[{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$. This is the approach adopted in ref. [1], and is also used here for the numerical computations.¹ In order to cover the whole single-particle spectrum, we have to choose $R\; \propto \; t$. In the limit $R\; \to \; \infty$, we recover the continuum case, but the discrete distribution also allows us to tackle special cases where only a few, or even just one frequency is better. For the remainder of the paper, we will focus on this discrete distribution.

As the cooling rate increases with $t$, longer cycle times lead to lower attainable energies. However, in the presence of noise, longer $t$ introduces more noise into the system. Therefore, for a given noise rate, there exists an optimal time $t$ that minimizes the final energy. The results of this optimization are summarized in table 1. We obtain that the optimal energy achievable under depolarizing noise scales as a power law with the noise rate $\kappa$. These results are derived analytically in the limits described above. We first derive the relative energy and then optimize it under those constraints (see table 1). This result shows that for the cooling algorithm, the energy scales with the noise rate as ${\kappa }^{1∕3}$ and does not depend on the quantum phase (see table 1). We also derive results on the achievable fidelity with this cooling process (see appendix A.1).

Adiabatic evolution. As is common in adiabatic algorithms, we consider a linear interpolation between the initial (simple) and the target Hamiltonian (see section 4.1). As before, increasing the total evolution time increases the performance of the algorithm. However, the effect of noise also increases with time. This interplay will dictate the optimal time and energy, leading to the results presented in table 1. Additionally, the performance of the adiabatic algorithm will strongly depend on whether the adiabatic path crosses the critical point at $𝜃\; =\; \pi ∕4$ or not. In case the target Hamiltonian still lies in the trivial phase (no crossing), the optimal time can be shown to scale like ${\kappa }^{-1∕3}$ and the achievable relative energy to scale like ${\kappa }^{2∕3}$ (see table 1). In the topological phase (), the gap is closed during the adiabatic evolution. Using Landau-Zener theory and including noise, we show that the relative energy scales as in eq. (5.6): we achieve for the optimal time a scaling proportional to ${\kappa }^{-2∕3}$ and the achievable energy as ${\kappa }^{1∕3}$ (see table 1).

It should be noted that scalings alone are not enough to determine the performance of the algorithm. The prefactors in front of them play an important role in the comparison between energy densities at finite values of $\kappa$. Taking them into account, it can be seen that the performance of the adiabatic algorithm in the selected range of $\kappa$ is better in the trivial phase than in the topological phase. That is, even though the rate at which relative energy grows with $\kappa$ is slower in the topological phase, the total relative energy is lower in the trivial phase. This can be understood from the absence of a level crossing in the trivial phase.

Numerical comparison between adiabatic and cooling algorithms. To complement the analytical results, we numerically optimize the runtime of the algorithms. This allows us not only to confirm the analytical predictions, but also to assess their performance beyond the assumptions of the analytical treatment. In figure 1 we present the difference between the numerically optimized (over the runtime) relative energies achieved with the adiabatic and the cooling algorithms as a function of the noise rate $\kappa$ and the Hamiltonian parameter $𝜃$. We perform this comparison for the initial state ${|0⟩}^{\otimes N}$ and the more realistic initial state ${\rho }^{cool}(\pi ∕2)$. As can be clearly seen, the adiabatic algorithm outperforms the cooling algorithm in the trivial phase. This is expected, as for $𝜃\; =\; \pi ∕2$ the ground state coincides with the initial state (or is close to it in the realistic scenario). However, cooling is more effective in the topological phase due to the gapless transition in the adiabatic path.

(a)

(b)

(c)

(d)

Figure 1: Difference $e_{opt}^{cool}-e_{opt}^{adiab}$ between the optimal relative energies achieved by the cooling and the adiabatic algorithms, plotted as a function of $𝜃∕\pi$ (horizontal axis) and ${\log <!--\; FUNCTION\; APPLICATION\; -->}_{10}\kappa$ (vertical axis) for $N\; =\; 100$ systems. Negative values (blue/dark) indicate that cooling gives lower energy (region labeled “C”); positive values (yellow/bright) indicate that the adiabatic algorithm gives lower energy (region labeled “A”). The black contour marks $e_{cool}^{opt}=e_{adiab}^{opt}$. (a) Topological phase, adiabatic algorithm starting from ${|0⟩}^{\otimes N}$ (ideal): cooling gives lower energy throughout most of this phase. (b) Trivial phase, ideal adiabatic algorithm: the adiabatic algorithm gives lower energy everywhere due to its small non-adiabatic prefactor and the ${\kappa }^{2∕3}$ power law in this gapped regime. (c) Topological phase, realistic adiabatic algorithm starting from ${\rho }^{cool}(\pi ∕2)$: the region where cooling yields lower energy expands. (d) Trivial phase, realistic adiabatic: the imperfect initial state reduces the adiabatic advantage, and cooling yields lower energy for $\kappa \; \gtrsim \; 10^{-3}$. For cooling the relative energies are computed numerically at $N\; =\; 100$, whereas for relative energies achieved by the adiabatic algorithm we use the large-$N$ analytic formulas to suppress the finite-size effect. The deviations at high $\kappa$ are due to the fact that we enter a regime where the here chosen discretization of $R$ matters (see appendix A.4).

2.3 Comparison between numerically optimized algorithms

In addition to cooling and adiabatic evolution, we consider QAOA as a third approach to ground-state preparation. To include it on equal footing in our comparison, we compare it to a numerically optimized version of our cooling algorithm (see [1]). Its performance is evaluated under the same noise model and both are optimized over the same number of parameters, ensuring a consistent and fair comparison. This allows us to assess how a finite-depth variational protocol competes with cooling across different phases.

QAOA. As is common in QAOA, we consider $p$ layers, each consisting of two unitaries generated by non-commuting Hamiltonians $H_M$ and $H_P$, which are both translationally invariant and quadratic. Consequently, the unitaries can again be decomposed into $2\; \times \; 2$ unitaries acting on each mode pair $(k,-k)$. Each unitary is then determined by a single parameter.

Comparison. First, it should be noted that in all regions considered here, the adiabatic algorithm performs worse than both the QAOA and the cooling algorithms. Therefore, we focus on comparing the latter two. In figure 2, we show the difference between QAOA and the cooling algorithm when both are optimized over 12 and 20 free parameters, respectively. For the cooling protocol, we choose a coupling range $r_c\in \; \{0.5,1.5\}$ [1] and allow two different cycles (each characterized by its own bath frequency $\mathrm{\Delta }$, cycle time $t$, and coupling parameters), giving a total of $8(r_c+\; 1)\; \in \; \{12,20\}$free parameters. For QAOA, we use $p\; =\; 4(r_c+\; 1)\; \in \; \{6,10\}$ layers to have the same number of parameters. As mentioned above, we choose as the initial state the state ${\rho }^{cool}(\pi ∕2)$. In case the parameters can be chosen exactly, QAOA typically achieves lower energies in the trivial phase ($𝜃\; >\; \pi ∕4$), while cooling gives lower energies in the topological phase (). This is expected, as preparing an SPT state from a trivial product state is generally difficult for QAOA. Neither algorithm dominates across the entire phase diagram.

Robustness. Since the optimal parameters cannot be implemented perfectly, we also analyze the robustness of the algorithms to small deviations in the parameters (see figure 2). As shown there, the cooling algorithm is more resilient to such fluctuations than the QAOA algorithm over a wide range of parameters.

2.4 Overall picture

Overall, our analytical and numerical results show that the optimal strategy is strongly phase-dependent: adiabatic evolution is advantageous in the trivial phase, whereas cooling is favorable in the topological phase. Our numerical results show that QAOA performs competitively with cooling in the trivial phase but is generally outperformed by cooling in the topological phase. We further find that cooling is more robust than QAOA to parameter imperfections, reinforcing its advantage in realistic settings.

3 Model

As in [1], we consider a family of translationally invariant quadratic fermionic Hamiltonians $H_S(𝜃)$ parametrized by $𝜃$. Both analytic calculations and large-$N$ numerics are feasible, and all dynamics decompose into independent momentum blocks. Below we define the Hamiltonian, the bath used for cooling, and the noise channel.

3.1 System

The system is a ring of $N$ free fermions (with $N$ even) with nearest-neighbor hopping and pairing, parametrized by $𝜃\; \in \; [0,\pi ∕2]$ [1]:

where $a_n$ ($a_n^{†}$) are fermionic annihilation (creation) operators with $\{a_n,a_m\}\; =\; \{a_n^{†},a_m^{†}\}\; =\; 0$, $\{a_n,a_m^{†}\}\; ={\delta }_{n,m}$, and we impose periodic boundary conditions $a_{N+1}=a_1$.

As mentioned before, this Hamiltonian has two phases separated by a critical point at $𝜃\; =\; \pi ∕4$. For $𝜃\; \to \; \pi ∕2$, the ground state is continuously connected to the trivial product vacuum (in the original fermionic basis). For $𝜃\; \to \; 0$, the model approaches the Kitaev/Majorana-chain limit (topological phase). At $𝜃\; =\; \pi ∕4$ the gap closes at momentum $k\; =\; N∕2$ and the ground state becomes degenerate. This mode is negligible for the energy density in the thermodynamic limit, but it becomes relevant for comparisons based on fidelity.

Since $H_S$ is quadratic and translation invariant, a Fourier transform decomposes it into independent $(k,-k)$ blocks. Each $(k,-k)$ block further splits into sectors of even and odd fermion parity: the even sector is spanned by $\{|{\mathrm{\Omega }}_S⟩,{ã}_k^{†}{ã}_{-k}^{†}|{\mathrm{\Omega }}_S⟩\}$ and the odd sector by $\{{ã}_k^{†}|{\mathrm{\Omega }}_S⟩,{ã}_{-k}^{†}|{\mathrm{\Omega }}_S⟩\}$, where $|{\mathrm{\Omega }}_S⟩$ is the particle vacuum of the system and ${ã}_k$ are the Fourier-transformed mode operators.

Because the ground state of $H_S(\pi ∕2)$ lies in the even-parity sector and unitary evolution generated by quadratic Hamiltonians preserves the fermion parity, both adiabatic evolution and QAOA reduce to independent effective qubits within this sector: the two-dimensional even-parity subspace of each $(k,-k)$ block defines a qubit with basis states $|0_k⟩\equiv |{\mathrm{\Omega }}_S⟩$ and $|1_k⟩\equiv {ã}_k^{†}{ã}_{-k}^{†}|{\mathrm{\Omega }}_S⟩$. For the cooling algorithm, which involves a CP map obtained by tracing out or resetting the bath after joint evolution with the system, the parity of only the system is not conserved in general (it is however conserved when considering system and bath together). Nevertheless, the evolution remains Gaussian and factorizes across momentum blocks, so that the analytic and numerical treatment still reduces to independent small subsystems of at most four dimensions. Throughout the remainder of this paper, when we refer to “mode $k$,” we mean the mode pair $(k,-k)$ and the corresponding effective qubit in the even-parity sector, unless stated otherwise.

Diagonalizing $H_S$ by a Bogoliubov transformation yields the mode Hamiltonians and energies

We call spectrum range the interval within which the energies of the modes vary, namely from ${𝜖}_m=\sqrt{1\; -\sin <!--\; FUNCTION\; APPLICATION\; -->2𝜃}$ at $k\; =\; N∕2$ to ${𝜖}_M=\sqrt{1\; +\sin <!--\; FUNCTION\; APPLICATION\; -->2𝜃}$ at $k\; =\; 0$. The gap closes (${𝜖}_m\to \; 0$) at the critical point $𝜃\; =\; \pi ∕4$.

As mentioned before, the 2-dimensional subspace of even parity is spanned by $\{|0_k⟩,|1_k⟩\}$. The Hamiltonian block for each mode reads then

$$H_k(𝜃)\; =w_k(𝜃){\sigma }_z+r_k(𝜃){\sigma }_x,$$

(3.4)

with $w_k(𝜃)\; =\sin <!--\; FUNCTION\; APPLICATION\; -->𝜃\; +\cos <!--\; FUNCTION\; APPLICATION\; -->𝜃\cos <!--\; FUNCTION\; APPLICATION\; -->{\varphi }_k$, $r_k(𝜃)\; =\cos <!--\; FUNCTION\; APPLICATION\; -->𝜃\sin <!--\; FUNCTION\; APPLICATION\; -->{\varphi }_k$, and ${\varphi }_k=\; 2\pi k∕N$. The total relative energy can be written in terms of the energies per mode as follows

$$e\; =\frac{\underset{k=0}{\overset{N∕2}{\sum }}{𝜖}_k{e}_k}{\underset{k=0}{\overset{N∕2}{\sum }}{𝜖}_k},$$

(3.5)

where $e_k=\; (E_k+{𝜖}_k)∕{𝜖}_k$ is the relative energy of mode $k$; here $E_k$ is the energy contribution of mode $k$ and $-{𝜖}_k$ is its ground-state value. For details (Bogoliubov angles, special cases $k\; =\; 0,N∕2$, etc.), we refer the reader to ref. [1].

3.2 Bath and Coupling

For the cooling algorithm, we introduce an auxiliary bath that can be prepared (and reset) close to its ground state. As in ref. [1], we consider $N$ independent fermionic modes with tunable onsite potential $\mathrm{\Delta }$:

$$H_B=\frac{\mathrm{\Delta }}{2}\underset{n=1}{\overset{N}{\sum }}\left(b_n^{†}{b}_n-b_n{b}_n^{†}\right),$$

(3.6)

where $b_n,b_n^{†}$ are fermionic operators acting on the bath.

The system modes ($a_n$) and bath modes ($b_n$) are coupled in the following translationally invariant way, up to a maximum range $r_c$:

$$V_{SB}=\; g\underset{n=1}{\overset{N}{\sum }}\sum _{j=-r_c}^{r_c}\left[\left({\lambda }_j{a}_n^{†}{b}_{n+j}+\; i{\mu }_j{a}_n{b}_{n+j}\right)+\text{h.c.}\right].$$

(3.7)

Here, $g$ is an overall coupling strength and ${\lambda }_j,{\mu }_j$ are dimensionless parameters.

In momentum space, after diagonalizing the Hamiltonian (in the Hilbert space corresponding to the system) via the Bogoliubov transformation, the joint Hamiltonian $H_{SB}=H_S+H_B+V_{SB}$ decomposes into independent $(k,-k)$ blocks. For generic $k$, one can write $H_{SB}={\sum }_k{\widehat{\alpha }}_k^{†}{h}_{SB,k}{\widehat{\alpha }}_k$ with ${\widehat{\alpha }}_k={({â}_k,{â}_{-k}^{†},{\widehat{b}}_k,{\widehat{b}}_{-k}^{†})}^{\top }$, and

where $A_k,B_k$ depend on the Bogoliubov angle ${\phi }_k$ and the couplings $\{{\lambda }_j,{\mu }_j\}$:

with ${\phi }_k=\frac{1}{2}{\tan <!--\; FUNCTION\; APPLICATION\; -->}^{-1}\left(\frac{\sin <!--\; FUNCTION\; APPLICATION\; -->\frac{2\pi k}{N}}{\tan <!--\; FUNCTION\; APPLICATION\; -->𝜃+\cos <!--\; FUNCTION\; APPLICATION\; -->\frac{2\pi k}{N}}\right)$.

3.3 Noise

We model noise by a depolarizing channel acting on both system and bath during each step of coherent evolution. In Lindblad form it is described by

where $\kappa$ denotes the noise strength. Since the Lindblad operators are site-local, this channel factorizes across sites, and also across modes; i.e. $e^{{\mathcal{\text{L}}}_{E}t}$ is the tensor product over all $e^{{\mathcal{\text{L}}}_E^{k}t}$, with the CPTP map [1]:

$$e^{{\mathcal{\text{L}}}_E^{k}t}({\rho }_k)\; ={\mathcal{\text{𝒟}}}_t({\rho }_k)\; =e^{-2\kappa t}{\rho }_k+\; (1\; -e^{-2\kappa t})\frac{𝟙}{4}.$$

(3.12)

This channel commutes with any unitary generated by a quadratic fermionic Hamiltonian (see [1]), which has distinct consequences for the three algorithms.

For the unitary algorithms (adiabatic and QAOA), the noise accumulated over the entire runtime can equivalently be collected into a single depolarizing channel ${\mathcal{\text{𝒟}}}_T$ applied after the complete unitary circuit; the energy of any quadratic observable is then damped by a factor $e^{-2\kappa T}$. For the cooling algorithm, the noise within each cycle can similarly be moved to act after that cycle’s unitary evolution step; however, the non-unitary bath reset between successive cycles prevents combining noise from different cycles into a single channel. As a result, increasing the total runtime of the algorithm degrades performance for all three algorithms, but the mechanism by which the optimal time balances coherent evolution against noise is algorithm-specific and constitutes the central analysis of section 5.

4 Quantum algorithms considered

Below we describe the three algorithms compared in this work: adiabatic evolution, Quantum Approximate Optimization Algorithm (QAOA), and cooling. For each algorithm, we also explain how depolarizing noise affects its performance.

4.1 Adiabatic algorithm

We consider the adiabatic preparation of the ground state of a target Hamiltonian $H_S({𝜃}_f)$ starting from the ground state of an initial Hamiltonian $H_S({𝜃}_i)$. We choose ${𝜃}_i=\; \pi ∕2$, corresponding to a trivial product state, and evolve towards a target $H_S({𝜃}_f)$. Following the standard adiabatic protocol [3, 35], we employ a linear interpolation between the initial and target Hamiltonians:

$$H(s)\; =\; (1\; -s)H_S({𝜃}_i)\; +\; s{H}_S({𝜃}_f),$$

(4.1)

where $s\; =\; t∕T\; \in \; [0,1]$ is the normalized time and $T$ is the total evolution time.

Due to the translational invariance the quadratic Hamiltonian $H(s)$ decomposes into independent $2\; \times \; 2$ blocks for each momentum mode $k$ (see appendix B for details). The many-body dynamics thus reduces to the independent evolution of $N∕2$ two-level systems. Depending on the targeted value of ${𝜃}_f$, the interpolation path may cross the quantum critical point at ${𝜃}_c=\; \pi ∕4$. We analyze how the relative energy scales in two regimes: (i) the gapped regime, described by adiabatic perturbation theory [3], where the deviation from the ground-state energy in the absence of noise scales as $T^{-2}$, and (ii) the gapless regime (crossing the critical point), described by the Kibble-Zurek mechanism [22, 23, 36–39] and Landau-Zener transitions [40, 41], where it scales as $T^{-1∕2}$.

In the ideal case, the initial state is the exact ground state ${|0⟩}^{\otimes N}$ at ${𝜃}_i=\; \pi ∕2$. In the realistic scenario, the initial state is the mixed state ${\rho }^{cool}(\pi ∕2)$ obtained from a cooling protocol at $𝜃\; =\; \pi ∕2$. Since ${\rho }^{cool}(\pi ∕2)$ is diagonal in the mode basis, the adiabatic evolution of each mode can be treated independently for both cases (see appendix B.4). Moreover, as the depolarizing channel [eq. (3.12)] commutes with the unitary evolution generated by $H(s)$, for every traceless quadratic observable $O$ and every state $\rho$ one obtains

$$tr\left[O{e}^{{\mathcal{\text{L}}}_{E}t}(\rho )\right]=e^{-2\kappa t}tr(O\rho ).$$

(4.2)

Therefore, the cumulative effect of noise over the total runtime $T$ is, at the level of quadratic observables such as the energy, equivalent to a single damping factor² applied after the ideal unitary evolution: the total energy of the target Hamiltonian is multiplied by $e^{-2\kappa T}$. The same argument applies to QAOA (see below). This creates a trade-off: increasing $T$ reduces deviations from the ground-state energy due to non-adiabatic transitions, but simultaneously enhances the effect of noise. The optimal runtime $T$, which balances these competing effects, is analyzed in section 5.

4.2 QAOA

The QAOA circuit consists of $p$ layers, each composed of two unitaries generated by the non-commuting Hamiltonians $H_M$ and $H_P$. We have

$$\mathcal{\text{𝒰}}(\mathbf{\text{𝛃}},\mathbf{\text{𝛄}})[\rho ]\; =\; U\rho U^{†},U\; =\prod _{l=1}^p{e}^{-i{\beta }_l{H}_M}{e}^{-i{\gamma }_l{H}_P},$$

(4.3)

where the parameters $\mathbf{\text{𝛃}},\mathbf{\text{𝛄}}\in {ℝ}^p$ are optimized to minimize the energy of $H_S({𝜃}_f)$. We take $H_M$ and $H_P$ to be translationally invariant and quadratic, so that the circuit factorizes into independent $2\; \times \; 2$ unitaries in each $(k,-k)$ block of even parity. Specifically, within the even-parity qubit basis defined in section 3.1, the single-mode generators read

$$h_{M,k}={\sigma }_z,h_{P,k}=\cos <!--\; FUNCTION\; APPLICATION\; -->{\varphi }_k{\sigma }_z+\sin <!--\; FUNCTION\; APPLICATION\; -->{\varphi }_k{\sigma }_x,$$

(4.4)

and the $l$-th QAOA layer acts on mode $k$ as

$$U_{k,l}({\beta }_l,{\gamma }_l)\; =e^{-i{\beta }_l{h}_{M,k}}{e}^{-i{\gamma }_l{h}_{P,k}}.$$

(4.5)

The full circuit factorizes as $U\; ={\otimes }_k{U}_k$ with $U_k={\prod }_{l=1}^p{U}_{k,l}$. While each $U_k$ is a product of $2\; \times \; 2$ unitaries, the parameters $\mathbf{\text{𝛃}},\mathbf{\text{𝛄}}$ are shared across all modes and must be optimized jointly to minimize the total energy.

Starting from an initial state ${\rho }_{init}={\otimes }_k{\rho }_{k,init}$, the variational energy is given by

$$E(\mathbf{\text{𝛃}},\mathbf{\text{𝛄}})\; =\underset{k=0}{\overset{N∕2}{\sum }}tr\left[U_k{\rho }_{k,init}{U}_k^{†}{H}_k({𝜃}_f)\right],$$

(4.6)

where $H_k({𝜃}_f)$ is the effective qubit Hamiltonian from eq. (3.4).

Using again that the depolarizing channel [eq. (3.12)] commutes with any unitary generated by a quadratic Hamiltonian, its effect on quadratic observables is equivalent to applying a single depolarizing channel ${\mathcal{\text{𝒟}}}_{t_{circ}}$ after the full circuit, where the total coherent time is $t_{circ}={\sum }_l(|{\beta }_l|\; +\; |{\gamma }_l|)$. Because each mode Hamiltonian $H_k$ is traceless, the energy of the maximally mixed state vanishes, and the energy of mode $k$ reduces to

$$E_{\kappa }=e^{-2\kappa t_{circ}}E(\mathbf{\text{𝛃}},\mathbf{\text{𝛄}}).$$

(4.7)

We use this expression for the numerical comparisons in section 6.

4.3 Cooling algorithm

A cooling algorithm extracts energy from the system by repeatedly coupling it to a bath prepared in its ground state, allowing joint evolution for a cycle time $\tau$, and then resetting the bath. In the absence of noise, each elementary cycle is described by the completely positive trace-preserving map

$$\mathcal{\text{E}}({\rho }_S)\; ={tr}_B\left[e^{-i{H}_{SB}\tau }({\rho }_S\otimes {\rho }_B)e^{i{H}_{SB}\tau }\right],$$

(4.8)

where ${\rho }_S$ is the system density matrix, ${\rho }_B$ denotes the initial state of the bath, prepared in its vacuum state, and ${tr}_B$ denotes the partial trace over the bath. The total Hamiltonian is $H_{SB}=H_S+H_B+V_{SB}$, as defined in section 3. More sophisticated schemes consist of sequences of multiple cycles with varying parameters, such as the system-bath couplings, the cycle time $\tau$, or the bath frequency $\mathrm{\Delta }$.

The analytical framework for our algorithm, developed in ref. [1], yields the following effective description. The analytical results are derived in the weak coupling limit, the low noise limit and the cooling limit (see appendix A for details). In the weak-coupling limit, ${(g\tau )}^2\ll \; 1$, a Dyson expansion of the map [eq. (4.8)] produces effective Lindblad dynamics for each momentum mode $k$. Two competing processes appear: a cooling process, peaked when the bath frequency $\mathrm{\Delta }$ matches the mode energy ${𝜖}_k$, and a heating process, which remains small in general provided the cooling limit ${𝜖}_k\tau \; \gg \; 1$ holds. A single bath frequency $\mathrm{\Delta }$ cools only modes near resonance, ${𝜖}_k\approx \mathrm{\Delta }$, and for a fixed cycle time, accidental resonances at $({𝜖}_k-\mathrm{\Delta })\tau \; =\; 2\pi n$ ($n\; \in \; \mathbb{Z}\; \setminus \; \{0\}$) lead to unwanted heating. This can be avoided using two different methods [1]: the first one randomizes the cycle times and adds multiple frequencies; the second couples each bath mode to more than one of the system’s modes and optimizes the cooling with respect to the corresponding couplings. The latter, in fact, does not need to operate in the weak coupling limit and gives rise to better performance for the models studied in [1], although one loses the analytical understanding of the first one. Below we detail how both methods operate.

Multi-time and frequency cooling. To cool all modes simultaneously, we employ a multifrequency strategy with local couplings (${\lambda }_0={\mu }_0=\; 1$ in eq. (3.7)). For each cooling cycle, we choose $\tau$ and $\mathrm{\Delta }$ randomly according to a given distribution, $D(\tau ,\mathrm{\Delta })\; =D_t^0(\tau )D_{R,{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M}^1(\mathrm{\Delta })$. The distribution is parametrized in terms of the mean cycle time, $t$, the number of frequencies, $R$, and the minimum and maximum values of the bath frequencies, ${\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M$, respectively, which depend on ${𝜖}_{m,M}$, the bandwidth of the system’s single-particle spectrum. Eventually, we also scale $R\; \propto \; t$, as explained below.

The distribution $D^0$ is uniform, namely $\tau \; \in \; [0,2t]$ with uniform probability distribution, so that $t$ is the averaged time of the cycle. The value of $t$ has to be chosen so that the weak coupling limit ${(gt)}^2\ll \; 1$ is fulfilled. The distribution $D^1$ is also uniform, but takes $R$ discrete values in the interval $[{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$, as explained previously.³ In order to choose those quantities, we note that in the weak coupling limit the cooling rate as a function of the mode energy ${𝜖}_k$ is given by a Lorentzian function centered at $\mathrm{\Delta }$ and with a linewidth ${\gamma }_0\propto \; 1∕t$ (see eq. (A.27) and figure 3). That is, if we choose one particular $\mathrm{\Delta }$ it will cool down the modes $k$ with ${𝜖}_k$ around the interval $\mathrm{\Delta }\pm {\gamma }_0$. Thus, we have to choose $R$ such that for all $k$, ${𝜖}_k$ lies within that interval for some ${\mathrm{\Delta }}_i$ (see figure 3 and appendix A.1.5 for details).

In order to derive analytical results, we average the map [eq. (4.8)] with respect to the distribution $D(\tau ,\mathrm{\Delta })$ and then expand the result in the weak coupling limit. In appendix A we show that after a sufficient number of cycles, the resulting state for each pair of modes $(k,-k)$ is $O(g^2{t}^2)$ close to the one obtained by those approximations. In practice, it is simpler to first expand in the weak coupling limit and then carry out the average, as it was done in [1].⁴ The result is that for each pair of modes, we obtain a rate equation with averaged cooling and heating rates, ${\gamma }_k^c$ and ${\gamma }_k^h$ (see appendix A). We also estimate, in appendix A.1.4, the total time needed to reach a target fidelity, finding an improvement over the bound of [1]. The relative energy of the steady-state of mode $k$ in the absence of noise is

$$e_k=\frac{2{\gamma }_k^h}{{\gamma }_k^c+{\gamma }_k^h}.$$

(4.9)

Depolarizing noise of strength $\kappa$ modifies this expression. Assuming a low noise limit ($\kappa t\; \ll \; 1$), each elementary subcycle of duration $\tau$ acquires a noise contribution $\kappa \tau$ [cf. eq. (3.12)]; averaging over the cycle time ($⟨\tau ⟩\; =\; t$) thus adds a term $\kappa t$ to both the cooling and heating rates, resulting in

$$e_k(t)\; =\frac{2{\gamma }_k^h+\; 2\kappa t}{{\gamma }_k^c+{\gamma }_k^h+\; 2\kappa t}.$$

(4.10)

When cooling dominates (${\tilde{\gamma }}_k^c\gg \; \kappa t$), the energy approaches the ground state ($e_k\to \; 0$); when noise dominates, the system depolarizes ($e_k\to \; 1$). The optimal trade-off between these regimes, and its scaling with $\kappa$, is derived in section 5.1. Let us note here that we could have also chosen a cooling algorithm where, at each step, $\tau$ and $\mathrm{\Delta }$ are chosen randomly and uniformly from the interval $[0,2t]$ and $[{𝜖}_m,{𝜖}_M]$ resp. (see appendix A). However, due to the fact that the discrete values might be experimentally more accessible, we proceed here with the discrete distribution of $\mathrm{\Delta }$, detailed in section 5.1.

Variationally optimized cooling. The second strategy we employ here uses a cooling procedure in which a cooling cycle is repeated until steady state is reached. Each global cycle is composed of a small number of cooling cycles with nonlocal couplings (coupling range $r_c$; see eq. (3.7)), each cycle having its own bath frequency $\mathrm{\Delta }$, cycle time $\tau$, and coupling parameters $\{{\lambda }_j,{\mu }_j\}$. All parameters are numerically optimized to minimize the total relative energy. As shown in ref. [1], this approach can achieve lower energies with shorter cycle times than the multifrequency strategy, at the cost of numerical fine-tuning. When comparing this algorithm to QAOA in section 6, we assess both its performance under noise and its sensitivity to imperfections in the optimized parameters.

5 Analytical treatment: multifrequency cooling vs. adiabatic algorithm

Having described the three algorithms, we now compare the performance of the multifrequency cooling and adiabatic approaches analytically. We determine the optimal relative energy each can achieve as a function of the noise rate $\kappa$. For the cooling algorithm, the analytical results are derived in three limits — the weak-coupling limit ${(gt)}^2\ll \; 1$, the cooling limit ${𝜖}_{m}t\; \gg \; 1$, and the low-noise limit $\kappa t\; \ll \; 1$, all of which are defined and used in appendix A; for these to be mutually compatible we require $\kappa \; \ll {𝜖}_m$.

In both cases, the optimal energy scales as a power law in $\kappa$, with an exponent that depends on the algorithm and the targeted quantum phase. Note that a lower exponent indicates greater resilience to noise since the energy degrades more slowly as $\kappa$ increases. We derive these exponents analytically for the multifrequency cooling and adiabatic algorithms below.

5.1 Multifrequency cooling: $e_{opt}^{cool}$ as a function of $\kappa$

The noisy steady-state energy given in eq. (4.10) depends on the cycle time $t$ through both the heating rates (which decrease with longer $t$ via sharper Lorentzian peaks) and the noise (which increases as $\kappa t$). In a noiseless setting, increasing $t$ improves spectral resolution, which causes $e_k$ to decrease monotonically with $t$. In a noisy setting a more careful optimization must be carried out, since the effect of noise will also increase with $t$. We want to optimize this tradeoff to determine the best achievable energy as a function of $\kappa$.

At fixed coupling strength $g$, the energy $e_k(t)$ in eq. (4.10) decreases monotonically with the cycle time $t$, but the limit $t\; \to \; \infty$is incompatible with the weak-coupling condition ${(gt)}^2\ll \; 1$ required for the perturbative treatment. We therefore impose that $gt$ is constant and set $gt\; =\; \eta$, fixing the product of coupling strength and cycle time. Numerical calculations confirm that the perturbative results remain accurate for $\eta \; \lesssim \; 0.5$; see appendix A and figure 9 therein. Under this constraint, the coupling strength decreases as $g\; =\; \eta ∕t$ and the number of bath frequencies grows as $R\; \simeq \sqrt{2∕3}({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)t$ (implemented as $R\; =\; \lfloor \sqrt{2∕3}({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)t\rfloor \; +\; 1$ when the integer nature of $R$ matters), leaving $t$ as the sole free parameter.

As shown in section 4.3, the averaged cooling rate ${\gamma }_k^c$ is $k$-independent, while the averaged heating rate ${\gamma }_k^h$ varies only weakly with $k$ and can be bounded from above (see appendix A and figure 7). Thus, the optimal cycle time is approximately the same for all modes. We can therefore optimize the total relative energy by setting $\partial e_k∕\partial t\; =\; 0$ at a representative mode and solving for the optimal time (see appendix A.2 for the full calculation). With the constraint that $gt$ is constant, in particular, $gt\; =\; \eta$, the noise-induced heating grows as $\kappa t$, while the averaged cooling rate decreases with growing $t$ (since ${\gamma }_k^c\propto g^{2}t\; =\; \eta ∕t$) and the heating rate decreases even faster with $t$ (since ${\gamma }_k^h\propto g^2={\eta }^2∕t^2$). The energy therefore first decreases and then increases with $t$, reaching a minimum at an intermediate optimal time $t_{opt}^{cool}\sim {(\kappa {\eta }^{-2})}^{-1∕3}$; substituting back into eq. (4.10) gives the minimum relative energy:

where ${\mathrm{\Delta }}_{M,m}$ are the largest and smallest bath frequencies, and $A$ is evaluated at $k\; =\; N∕4$, i.e. in the middle of the spectrum (the weak $k$-dependence is discussed in appendix A).

5.2 Adiabatic algorithm: $e_{opt}^{adiab}$ as a function of $\kappa$

To compare the performance of the cooling algorithm with the adiabatic algorithm, we now derive the relative energy as a function of the noise rate for the latter. Starting from ${𝜃}_i=\; \pi ∕2$ and targeting ${𝜃}_f$, two qualitatively different regimes arise depending on whether the interpolation path [eq. (4.1)] crosses the quantum critical point at ${𝜃}_c=\; \pi ∕4$ or not. Since $H(s)$ is quadratic and translationally invariant, the many-body dynamics factorizes into independent two-level problems, one for each momentum mode $k$.

While each individual mode undergoes standard two-level dynamics, the nontrivial step is summing over all $N∕2$ modes — each with a different instantaneous gap and transition amplitude — to obtain the total energy and its scaling with $T$. Recall that, as the depolarizing channel commutes with the unitary evolution (see section 3.3), its effect on the energy is given by an overall damping factor $e^{-2\kappa T}$. Non-adiabatic effects increase the deviation from the ground-state energy for faster processes, whereas depolarizing noise leads to higher energies for longer evolution times. The competition between these two mechanisms gives rise to an optimal runtime that minimizes the achievable energy.

We sketch the derivations below (full details are given in appendix B). Throughout, we use the superscripts “prod” (product/trivial phase) and “ent” (entangled/topological phase) to distinguish the results for the two phases, as they give different scalings. In table 1 we summarize the scalings of the relative energy as a function of $\kappa$ for both the cooling and adiabatic algorithms.

Trivial phase, no crossing (${𝜃}_f>\; \pi ∕4$). When the adiabatic path does not cross the critical point, the spectral gap remains of order $\mathcal{\text{𝒪}}(1)$ for all modes throughout the evolution. Applying first-order adiabatic perturbation theory [42] to each $2\; \times \; 2$ mode block, we expand the time-evolved state in the instantaneous eigenbasis and find that the transition amplitude to the excited state is given by an oscillatory integral whose rapid phase oscillations (at rate $\propto \; T$) suppress contributions from the bulk of the integration domain; only the boundary terms survive. This yields transition probability $|c_1(k){|}^2\propto T^{-2}$ for each mode (see appendix B). Summing over all modes, the final energy in the absence of noise takes the form

$$E\; \approx E_{GS}\left(1\; -\frac{Ã}{T^2}\right),\text{ with }Ã\; \approx \frac{{\cos <!--\; FUNCTION\; APPLICATION\; -->}^2{𝜃}_f}{8}.$$

(5.2)

Including the damping factor due to noise, $e^{-2\kappa T}$, and expanding for $\kappa T\; \ll \; 1$, the total energy becomes $E\; \approx E_{GS}(1\; -Ã∕T^2-\; 2\kappa T)$. The non-adiabatic contribution ($\propto T^{-2}$) is suppressed at large $T$, while the noise contribution ($\propto \; \kappa T$) is small at short $T$. Optimizing over $T$ under the constraint $\kappa T\; \ll \; 1$ gives

$$T_{opt}^{prod}\approx {\left(\frac{Ã}{\kappa }\right)}^{1∕3},\text{ and }{e}_{opt}^{prod}\approx \; 3{Ã}^{1∕3}{\kappa }^{2∕3}.$$

(5.3)

Topological phase, crossing the critical point (). When the adiabatic path crosses the critical point, the spectrum of mode $k\; =\; N∕2$ becomes gapless, while nearby modes also develop very small gaps. As a result, the system is particularly susceptible to excitations in this region.

Close to the critical point, the dynamics of modes near $k\; =\; N∕2$ can be approximated by the well-known Landau-Zener model [40, 41]. This model calculates the excitation probability of a two-level system evolution that starts in the ground state and goes through an avoided crossing at some speed (for details, see appendix B). We use this model to analyze each mode $k$ in our system. The energy levels approach the avoided crossing with a “velocity” set by $v_{LZ}=\; 2∕T$, while the minimum gap in this crossing depends on how far the mode is from $k\; =\; N∕2$, which is captured by the parameter ${\delta }_k=\; 2\pi k∕N\; -\pi$. Within this approximation, the excitation probability takes a Gaussian form,

$$P_{LZ}^k=\exp <!--\; FUNCTION\; APPLICATION\; -->\left(-\pi GT{\delta }_k^2\right),$$

(5.4)

where $G\; ={\cos <!--\; FUNCTION\; APPLICATION\; -->}^2{𝜃}_f∕{(1\; +\cos <!--\; FUNCTION\; APPLICATION\; -->{𝜃}_f-\sin <!--\; FUNCTION\; APPLICATION\; -->{𝜃}_f)}^3$ depends only on the final parameter ${𝜃}_f$. For large evolution time $T$, this Gaussian becomes sharply peaked around $k\; =\; N∕2$, meaning that only a narrow range of modes is excited. The width of this region scales as $\sim \; 1∕\sqrt{T}$, which is the characteristic behavior predicted by the Kibble-Zurek mechanism [22, 23, 38].

By approximating the sum over modes as an integral, one finds that the energy deviation from the ground state decreases as

$$E\; \approx E_{GS}\left(1\; -\frac{\tilde{B}}{\sqrt{T}}\right),$$

(5.5)

where $\tilde{B}$ is a constant of order one that depends on ${𝜃}_f$ and on the behavior of the spectrum near the critical point.

Including the effect of noise and optimizing over the total evolution time $T$ yields

$$T_{opt}^{ent}\approx {\left(\frac{\tilde{B}}{4\kappa }\right)}^{2∕3},e_{\text{opt}}^{ent}\approx \frac{3}{2^{1∕3}}{\tilde{B}}^{2∕3}{\kappa }^{1∕3}.$$

(5.6)

Compared to the trivial phase, the energy approaches the ground state more slowly with increasing $T$ (as $T^{-1∕2}$ rather than $T^{-2}$). Even though we find a more favorable scaling with $\kappa$ than in the trivial phase, the relative energy is higher due to the prefactors.

Finally, at very low noise, the energy cannot be reduced arbitrarily. This is because the mode $k\; =\; N∕2$ has exactly zero gap at the critical point and is therefore always excited, no matter how slow the evolution is. This sets a lower bound on the energy that is independent of $\kappa$. The crossover between this finite-size limit and the ${\kappa }^{1∕3}$ scaling occurs when ${\kappa }^{1∕3}\sim {𝜖}_{N∕2}∕|E_{GS}|$, and shifts to smaller values of $\kappa$ as the system size $N$ increases (see appendix B and figure 5).

Figure 4 compares the analytic cooling and adiabatic algorithms across four representative target values of ${𝜃}_f$, all starting from ${𝜃}_i=\; \pi ∕2$ on a system of $N\; =\; 100$ sites. For each value of $\kappa$, both the total adiabatic evolution time $T$ and the cooling cycle time $t$ are independently optimized to minimize the relative energy. One of the key findings is that the results are phase-dependent. In the trivial phase [panel (a), ${𝜃}_f=\; \pi ∕3$], the adiabatic ${\kappa }^{2∕3}$ scaling gives lower energy than cooling’s ${\kappa }^{1∕3}$ across the entire range of $\kappa$ shown, because the small prefactor $Ã\; \approx {\cos <!--\; FUNCTION\; APPLICATION\; -->}^2{𝜃}_f∕8$ keeps the adiabatic energy well below the cooling curve. In the topological phase [panels (b)–(d)], Landau-Zener excitations near the gap-closing momentum reduce the noiseless convergence to $T^{-1∕2}$, so after optimization both algorithms scale as $e\; \propto {\kappa }^{1∕3}$. The relative performance is then determined by prefactors: as ${𝜃}_f$ decreases deeper into the topological phase, the adiabatic prefactor $\tilde{B}$ grows while the cooling prefactor stays roughly constant. The ideal adiabatic algorithm gives slightly lower energy at ${𝜃}_f=\; \pi ∕6$ [panel (b)], the two are comparable at ${𝜃}_f=\; \pi ∕8$ [panel (c)], and cooling gives lower energy at ${𝜃}_f=\; \pi ∕10$ [panel (d)].

At small $\kappa$ in panels (b)–(d), the adiabatic curves flatten to a noise-independent floor: the mode $k\; =\; N∕2$ has exactly zero gap at the critical point and is unavoidably excited regardless of sweep time. This floor shifts to smaller $\kappa$ with increasing $N$ (figure 5).

Figure 1 maps out the difference $e_{cool}-e_{adiab}$ across the full $(𝜃,\kappa )$ plane in the limit $N\; \gg \; 1$. Dark regions indicate that cooling achieves lower energy (“C”); bright regions (“A”) indicate that the adiabatic algorithm does; the black contour marks equal performance. In the topological phase [panels (a, c)], cooling outperforms the adiabatic algorithm throughout most of the plane. In the trivial phase [panel (b)], the adiabatic algorithm, with scaling ${\kappa }^{2∕3}$ leads to lower energies. Panels (c, d) show the realistic scenario (initial state ${\rho }^{cool}(\pi ∕2)$); the imperfect initialization adds a ${\kappa }^{2∕3}$ contribution that shifts the equal-performance contour in favor of cooling, especially at high noise [panel (d), $\kappa \; \gtrsim \; 10^{-3}$].

$\pi ∕2\; \to \pi ∕3$

(a)

$\pi ∕2\; \to \pi ∕6$

(b)

$\pi ∕2\; \to \pi ∕8$

(c)

$\pi ∕2\; \to \pi ∕10$

(d)

Figure 4: Optimal relative energy $e_{\text{opt}}$ (log scale) as a function of the noise strength $\kappa$ (log scale) for $N\; =\; 100$ sites, comparing three protocols: analytic cooling with $gt\; =\; 0.5$ and ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}\pm \; 0.4$ (blue circles), ideal adiabatic with initial state ${|0⟩}^{\otimes N}$ (orange squares), and realistic adiabatic with initial state ${\rho }^{cool}(\pi ∕2)$ (green squares). Dashed lines show the analytic predictions from eq. (5.1) (blue, ${\kappa }^{1∕3}$) and appendix B (orange, green). (a) ${𝜃}_f=\; \pi ∕3$ (trivial phase): the adiabatic ${\kappa }^{2∕3}$ scaling gives lower energy than the achievable rate of the cooling algorithm, which is ${\kappa }^{1∕3}$. (b) ${𝜃}_f=\; \pi ∕6$ (topological phase): both scale as ${\kappa }^{1∕3}$; the adiabatic algorithm gives slightly lower energy. (c) ${𝜃}_f=\; \pi ∕8$: the two algorithms become comparable. (d) ${𝜃}_f=\; \pi ∕10$: cooling gives lower energy than both adiabatic protocols. In panels (b)–(d), the adiabatic curves flatten at small $\kappa$ due to the unavoidable excitation of the zero-gap mode $k\; =\; N∕2$.

These power-law behaviors found analytically rely on the assumption of large-$N$; for sufficiently small $\kappa$, the adiabatic curves can cross over to a different regime due to finite-size effects (see figure 5), so the ${\kappa }^{2∕3}$ and ${\kappa }^{1∕3}$ laws found analytically hold within a finite window $\kappa \; \in \; [{\kappa }_{min},{\kappa }_{max}]$, as shown in appendix B.

6 Numerical optimization: QAOA vs optimized cooling

The analytical expressions for the relative energy as a function of the noise rate, derived in section 5, enable a fair comparison between the multifrequency cooling and adiabatic algorithms. Here, we study how the variational approach, QAOA, performs relative to those state-preparation protocols. We perform a direct numerical optimization over the full parameter space $(𝜃,\kappa )$.

For the variational protocols studied in section 6, i.e. QAOA and optimized cooling, we find that both outperform the particular linear-ramp adiabatic protocol of eq. (4.1). Although alternative adiabatic paths could yield improved performance, optimizing such paths is beyond the scope of this work. Consequently, the adiabatic algorithm is excluded from the final comparison performed here.

Our numerical comparison proceeds in two steps. First, we compare the relative energies achieved by cooling and QAOA using optimally tuned parameters. Since such precise parameter settings are typically not accessible in experiments, we then analyze the robustness of both approaches under perturbations of the optimized parameters.

Before we begin with the comparison, let us make the following general remark. While QAOA performs well in the trivial phase, it struggles to prepare the ground state near $𝜃\; =\; 0$. The limit $𝜃\; \to \; 0$ corresponds to a Kitaev/Majorana chain, which realizes a symmetry-protected topological (SPT) phase [43]. Preparing an SPT state from a trivial product state using a local unitary circuit such as in QAOA requires a circuit depth that scales linearly with $N$ [44, 45]. Since the QAOA circuits considered here have constant depth in the sense that the number of layers $p$ does not scale with $N$, and the initial state is close to a product state, QAOA is expected to face difficulties in efficiently preparing a state with high fidelity with respect to the ground state. This result gives some intuition as to what energies we should expect to observe in both phases for QAOA. In contrast to that, the cooling algorithm is non-unitary and thus not subject to these circuit-depth constraints (other cooling and dissipative approaches based on continuous sweeps have also been explored for topological phases [49, 50]). Since the performance of the cooling algorithm mostly depends on the energy spectrum of the system, and this spectrum is symmetric with respect to $𝜃\; =\; \pi ∕4$, we expect cooling to work roughly equally well in both phases.

6.1 Optimal energies with exact optimal parameters

This subsection compares the best energies obtained when QAOA and optimized cooling are implemented with their numerically optimized parameters.

We proceed with the comparison as follows. For the cooling protocol, we choose a coupling range $r_c$ and allow two different cooling cycles. Each of them is characterized by its own bath frequency $\mathrm{\Delta }$, cycle time $t$, and coupling parameters $\{{\lambda }_j,{\mu }_j\}$, giving a total of $N_{pars}=\; 8(r_c+\; 1)$ free parameters. For QAOA, we use $p\; =\; 4(r_c+\; 1)$ layers, giving the same total number of free parameters, $2p$, as we have two parameters ${\beta }_l,{\gamma }_l$ for each layer. For the numerics, we chose $r_c=\; 0.5$ and $r_c=\; 1.5$, which amounts to 12 and 20 free parameters, respectively. The half-integer coupling ranges follow the convention of ref. [1], where $r_c=\; 0.5$ means only the $j\; =\; +1$ neighbor is coupled and not the $j\; =\; -1$.

Since the exact ground state ${|0⟩}^{\otimes N}$ at $𝜃\; =\; \pi ∕2$ cannot be prepared perfectly in the presence of noise, we use the same realistic initial state for both algorithms: the steady state obtained from cooling at $𝜃\; =\; \pi ∕2$, denoted ${\rho }^{cool}(\pi ∕2)$. Note that this steady state may, in general, differ from the one derived in appendix A.4 and utilized before, since the cooling protocol used here is the variationally optimized one rather than the analytic single-frequency scheme. To avoid inconsistencies, we compute this “realistic initial state” numerically for each parameter set.

All parameters are optimized using the BFGS algorithm as implemented in the Optim.jl package for Julia. For each $(𝜃,\kappa )$ grid point, we perform three independent optimization runs from random initial parameter vectors and retain the set yielding the lowest energy. We note that while we focus here on standard implementations with equal parameter counts, QAOA and cooling protocols could both be further improved using advanced optimization strategies [51, 52].

Figure 6 shows the optimal relative energies as three-dimensional surface plots over the $(𝜃∕\pi ,{\log <!--\; FUNCTION\; APPLICATION\; -->}_{10}\kappa )$ plane for $N\; =\; 20$ sites. Both surfaces decrease (i.e., improve) as $\kappa$ decreases and as $𝜃$ approaches $\pi ∕2$. The main difference arises near $𝜃\; =\; 0$ (topological phase). There, QAOA faces limitations, which might be related to the SPT depth constraint mentioned above, whereas cooling attains lower energies. In the trivial phase ($𝜃\; \to \; \pi ∕2$), the performance of both algorithms is comparable. Therefore, in the SPT phase, there is a clear advantage of cooling with respect to QAOA.

6.2 Parameter resilience

In practice, control parameters cannot be implemented exactly. To quantify the sensitivity of each algorithm to such imperfections, we introduce random relative perturbations to the optimized parameters. Specifically, for an optimized parameter vector $\mathbf{\text{𝜗}}=\; ({𝜗}_1,\dots ,{𝜗}_M)$, we sample i.i.d. variables $u_j\sim \mathrm{\text{Unif}}[-\delta ,\delta ]$ and replace each parameter as ${𝜗}_j\mapsto {𝜗}_j(1\; +u_j)$. For each perturbation strength $\delta$, we compute the resulting relative energy and average over 100 independent realizations.

Figure 2 shows the difference of the relative energies, $e_{\text{opt}}^{\text{cool}}-e_{\text{opt}}^{\text{QAOA}}$, across the $(𝜃,\kappa )$ plane. Panels (a, d) show the comparison with exact optimal parameters ($\delta \; =\; 0$): QAOA gives lower energies in the trivial phase ($𝜃\; >\; \pi ∕4$), while cooling gives lower energy in the topological phase (). Panels (b, e) show the same comparison with $2\%$ relative deviations ($\delta \; =\; 0.02$), and panels (c, f) with $5\%$ deviations ($\delta \; =\; 0.05$). As the deviations grow, the region in which cooling gives lower energy expands noticeably, demonstrating that the cooling protocol is more robust to parameter imperfections than QAOA.

This greater robustness of cooling can be attributed to its convergence to a fixed point: the steady state is an attractor of the cooling map, so small perturbations in the coupling parameters shift the attractor only slightly. In contrast, QAOA produces a single time-evolved state whose energy is directly sensitive to each parameter value, which does not have any stabilizing effect from a fixed-point structure.

7 Conclusion and outlook

We have developed a unified framework for comparing ground-state preparation algorithms — multifrequency cooling, adiabatic evolution, and QAOA — in the presence of depolarizing noise, applied to an exactly solvable family of quadratic fermionic Hamiltonians.

Our central analytical result is the derivation of scaling laws for the achievable relative energy as a function of the noise rate $\kappa$ (table 1). In the trivial (gapped) phase, the adiabatic algorithm achieves a ${\kappa }^{2∕3}$ scaling which, together with its small prefactor in the gapped regime, yields lower relative energies than the ${\kappa }^{1∕3}$ cooling result in the parameter range studied here. In the topological phase, where the adiabatic path crosses a quantum critical point, both algorithms scale as ${\kappa }^{1∕3}$. However, cooling achieves lower prefactors for the targeted values of ${𝜃}_f$ deep in the topological phase. No single algorithm dominates across the full phase diagram. That is, the optimal strategy is phase dependent.

Numerical optimization of QAOA and numerically optimized cooling algorithms lead to a similar picture. QAOA performs competitively with cooling in the trivial phase but is outperformed in the topological regime. This is consistent with the circuit-depth constraints on preparing symmetry-protected topological states from product states. The cooling protocol also exhibits greater robustness to parameter imperfections (figure 2), a practical advantage for realistic implementations.

Several open questions remain. Most importantly, future work should investigate whether the phase-dependent advantage of cooling persists for interacting, non-integrable Hamiltonians, where the structure of the spectrum and the nature of the phase transition change qualitatively. Optimizing the adiabatic schedule beyond the linear ramp (e.g., via counterdiabatic driving [46]) and combining cooling with variational optimization of bath parameters, or using machine learning methods to dynamically optimize the cooling cycle, are additional directions worth exploring.

Let us also mention that our multifrequency cooling strategy shares conceptual similarities with the filtered Lindbladian approach of Zhan et al. [12], though the mechanism (cooling vs dissipative state preparation), and also the implementations differ: static couplings with randomized cycle times and frequencies (ours) versus time-dependent couplings (theirs). A systematic comparison of these complementary approaches in terms of convergence time and time to solution would be valuable, especially in a noisy setting.

On the experimental side, incorporating Trotterization errors [47] for digital implementations and gate-level noise models would connect our continuous-time analysis to near-term superconducting [14] and trapped-ion [48] platforms. More broadly, the analytical approach developed here — deriving closed-form noise scaling laws and comparing them across quantum phases — provides a template for benchmarking ground-state preparation algorithms beyond the specific model and noise channel studied in this work.

Acknowledgments

This research is part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. The research is partly funded by THEQUCO as part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. The work is partially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2111 - 390814868. The work at MPQ is supported from the German Federal Ministry of Education and Research (BMBF) through the funded project ALMANAQC, grant number 13N17236 within the research program “Quantum Systems”. Furthermore, we acknowledge the BMW endowment fund. D.M. and B.K. acknowledge funding from the HORIZON-CL4-2021-DIGITAL-EMERGING-02-10 under grant agreement No. 101080085 (QCFD).

Appendix

A Multifrequency cooling: analytical derivations

In this appendix, we derive the analytical scaling laws for the multifrequency cooling protocol under depolarizing noise. We proceed in three stages: we first review the multifrequency cooling algorithms in the weak coupling regime, where we focus on the specific choice of times and frequencies presented in section 5 (appendix A.1); we then derive the steady-state relative energy of each mode in the presence of noise (appendix A.2) and optimize the cycle time to obtain the lowest achievable energy as a function of the noise rate (appendix A.3). Finally, in appendix A.4, we treat a single frequency protocol and show that it is advantageous in the edges of the spectrum. We also bridge the region close to the edges of the spectrum with a few-frequency protocol.

A.1 Weak-coupling and cooling limits

In section 5 we provided some analytical expressions for the energy in steady state after the cooling procedure in the presence of noise. Those are valid in the weak coupling and cooling limits, and for low noise rate. Furthermore, they will be obtained after averaging over different realizations of the cooling maps. In this subsection we show that such average faithfully describes the cooling process, and derive the analytical expressions for the steady state in both the weak coupling and cooling limits. In the following subsection, we will also include noise in this derivation.

We consider the map ${\mathcal{\text{E}}}_{k,\tau ,\mathrm{\Delta }}$ describing a cooling cycle for the density operator in a sector including the modes labeled by $(k,-k)$,⁵ where $\tau$ and $\mathrm{\Delta }$ are the cycle time and bath frequency, respectively. They are randomly chosen from flat distributions $\tau \; \in \; [0,2t]$ and $\mathrm{\Delta }\in \; [{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$ with

$${\mathrm{\Delta }}_M>{𝜖}_M\ge {𝜖}_m>{\mathrm{\Delta }}_m>\; 0,$$

where ${𝜖}_{m,M}$ are the minimum and maximum values of the mode energy $𝜖$. This ensures that all modes are cooled, as described in [1]. Note that we have assumed that the Hamiltonian is gapped, i.e. ${𝜖}_m>\; 0$.

For a given realization we can write the action of the next cooling cycles as a sequence

$${\rho }_{\ell +1}=E_{k,{\tau }_{\ell +1},{\mathrm{\Delta }}_{\ell +1}}({\rho }_{\ell })$$

(A.1)

where ${\rho }_{\ell }$ is the density operator corresponding to the modes with momenta $\pm \; k$ at cycle $\ell$, and ${\tau }_{\ell },{\mathrm{\Delta }}_{\ell }$ the corresponding values chosen for the map. In the following, for the sake of simplicity we will omit all the subscripts whenever they are not needed. We are interested in the state ${\rho }_{\ell }$ for $\ell \; \gg \; 1$. Although no steady state will be reached by applying all those maps, since at each time step we apply a different random channel, we will demonstrate that the final state (for $\ell \; \gg \; 1$) will be close to the steady state of the map that is obtained by averaging $E_{k,\tau ,\mathrm{\Delta }}$ with respect to $\tau$ and $\mathrm{\Delta }$, as long as we restrict ourselves to the weak-coupling limit $gt\; \ll \; 1$. Furthermore, we will use the cooling limit ${𝜖}_{m}t\; \gg \; 1$ in order to obtain simple analytical expressions for that steady state.

We have divided the remainder of this subappendix into five parts. In the first, we set up the method used to compute both the steady state of the averaged map and its standard deviation. From there it will be clear that averaging the map and expanding it in the weak coupling and cooling limits can be performed in any order to approximate the steady state of such a map. In the second, we explicitly compute that steady state in the weak-coupling limit. In the third, we evaluate the standard deviation and show that, in this limit, each individual realization — the state obtained after applying $\ell \; \gg \; 1$ cooling maps with random $\tau$ and $\mathrm{\Delta }$ — is close to the stationary state of the averaged map. In section A.1.4 we further expand the resulting state, and the cooling and heating rates, in the cooling limit. Finally, in section A.1.5 we consider a discrete distribution for the values of $\mathrm{\Delta }$, which is needed to analyze the dependence of the final energy on the number of frequencies, as discussed in section 5.1.

A.1.1 Averaged map and steady state

Given the fermion parity superselection rule, all the information about $\rho$ is contained in the following expectation values:

We can define a vector $P\; =\; (n,d,p,p^{\ast },q,q^{\ast })$, so that the sequence in eq. (A.1) can be equally written as

$$P_{\ell +1}=X_{\ell +1}{P}_{\ell }+Y_{\ell +1}$$

(A.3)

where $X_{\ell },Y_{\ell }$ are a matrix and a vector, respectively, which depend on both ${\tau }_{\ell }$ and ${\mathrm{\Delta }}_{\ell }$.

We are interested in $P_{\ell }$ in the limit $\ell \; \to \; \infty$, since this contains the information about the final energy and therefore about the success of the cooling process. For each realization, we will obtain a different result since it contains different choices of the cycle times and bath frequencies, and thus $P_{\ell }$ follows a distribution itself. In order to characterize it, we may consider its mean value and the standard deviation (or, more generally, the covariance matrix). In case the latter is small, the cooling will be well characterized by the first. In the next subsection we will show that the standard deviation is $O({(gt)}^2)$, so that in the weak-coupling limit we can confidently extract the conclusions about the success of the cooling process from the averaged values.

We define ${\overline{P}}_{\ell }$ as the vector obtained when averaging over $P_{\ell }$ with respect to all previous choices of $\tau$ and $\mathrm{\Delta }$, following the distributions mentioned above. Then we have,

$${\overline{P}}_{\ell +1}=\overline{X}{\overline{P}}_{\ell }+\; Ȳ$$

(A.4)

where the bar indicates the average with respect to the $\ell \; +\; 1$ random values of $\tau$ and $\mathrm{\Delta }$ chosen in the $\ell \; +\; 1$ different steps. Note that the average of $X_{\ell +1}{P}_{\ell }$ can be written as the product of averages since the first only depends on the choice of $\tau$ and $\mathrm{\Delta }$ in step $\ell \; +\; 1$, while the second depends on all previous steps. In the limit $\ell \; \to \; \infty$, we obtain

$$\overline{P}=\underset{\ell \to \infty }{lim}{\overline{P}}_{\ell }={(𝟙\; -\overline{X})}^{-1}Ȳ$$

(A.5)

In order to compute the standard deviation of the distribution, we consider

$$P_{\ell +1}\otimes P_{\ell +1}=\; (X_{\ell +1}{P}_{\ell }+Y_{\ell +1})\; \otimes \; (X_{\ell +1}{P}_{\ell }+Y_{\ell +1})$$

(A.6)

Then, we average with respect to all realizations, and take the limit $\ell \; \to \; \infty$ to obtain

$$Q\; =\overline{P\; \otimes P}-\overline{P}\otimes \overline{P}$$

(A.7)

which contains all the information about the corresponding covariance matrix. Once we obtain the results for $\overline{P}$ and $Q$, we can perform the expansions in the weak coupling and cooling limit, which involve expanding $\overline{X}$, $\overline{X\; \otimes X}$, etc. in terms of $gt$ and ${𝜖}_{m}t$. Crucially it holds that, for all those quantities, the average of the expansion coincides with the expansion of the average. Thus, we can first do the expansion in $gt$, then take the average, and finally the expansion in ${𝜖}_{m}t$. Indeed, we will do that in the subsequent subsections since this heavily simplifies the derivations.

A.1.2 Average in the weak-coupling limit

We will first expand the map $E_{k,\tau ,\mathrm{\Delta }}$ up to second order in $\eta \; =\; gt$ and then average it with respect to $\tau$ and $\mathrm{\Delta }$. Both steps were carried out in [1], so we will take some of the results from there. Note, however, that in that reference we averaged over a discrete distribution in $\mathrm{\Delta }$, while here we will consider a continuous distribution in order to get simpler formulas. At the end of this section we will also consider a discrete distribution in which the discretization depends on $t$, since most of the computations in the main text have been carried out for such a distribution.

In the weak coupling limit, ${\eta }^2={(gt)}^2\ll \; 1$ we have that the cooling map can be transformed via a Dyson expansion into the following Lindbladian form [1]:

Here, ${\mathcal{\text{L}}}_O$ denotes the Lindblad superoperator ${\mathcal{\text{L}}}_O(\rho )\; =\; O\rho O^{†}-\frac{1}{2}\{O^{†}O,\rho \}$, and

with $𝜖$ the energy of the fermionic modes and coefficients $x$ and $y$ given by

We have that $|x|,|y|\; \le \; 1$, since $\tau \; \in \; [0,2t]$ and the bound is reached for the resonances $(\mathrm{\Delta }\pm 𝜖)\tau \; \ll \; 1$. As explained in [1], the interpretation of this map is that apart from the system evolution, $e^{-i{h}_S\tau }$, we have dissipative terms resulting from the weak coupling to the bath.

We can now use this equation to obtain the change of the parameters in eq. (A.2) in one cooling cycle,

where we have defined $|z{|}^2=\; |x{|}^2+\; |y{|}^2$, and we have dropped the subscript $\ell$ in $x$, $y$ and $z$. From these equations, we can easily retrieve the matrix $X$ and the vector $Y$ defined in eq. (A.3) and then compute the corresponding averages in the limit $\ell \; \to \; \infty$. However, before doing that, note that for any realization we have that $|d_{\ell }|,|p_{\ell }|,|q_{\ell }|\; \to \; 0$ in that limit.⁶ Since we are neglecting terms $O({\eta }^4)$, we conclude that all those quantities are of that order in the limit $\ell \; \to \; \infty$. In particular, we can replace those values of $p_{\ell }$ in the equation for $n_{\ell +1}$ so that for sufficiently large $\ell$ we have

$$n_{\ell +1}=n_{\ell }+{\eta }^2\left[|y{|}^2-|z{|}^2{n}_{\ell }\right]+O({\eta }^4).$$

(A.15)

We can now average eq. (A.15) with respect to $\tau$ and $\mathrm{\Delta }$, obtaining

$${\overline{n}}_{\ell +1}-{\overline{n}}_{\ell }=\; -({\gamma }^c+{\gamma }^h){\overline{n}}_{\ell }+{\gamma }^h+O({\eta }^4),$$

(A.16)

where

$${\gamma }^{c,h}=\frac{g^2}{{\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m}\left[F({\mathrm{\Delta }}_M\mp 𝜖,t)\; -F({\mathrm{\Delta }}_m\mp 𝜖,t)\right],$$

(A.17)

with

$$F(𝜖,t)\; =\frac{\sin <!--\; FUNCTION\; APPLICATION\; -->(2𝜖t)}{2{𝜖}^{2}t}+\frac{\cos <!--\; FUNCTION\; APPLICATION\; -->\left(2𝜖t\right)-\; 2}{𝜖}+\; 2t\mathrm{\text{Si}}\left(2𝜖t\right)$$

and $\mathrm{\text{Si}}(y)$ is the sine integral function. We can thus write

$${\overline{n}}_k=\underset{\ell \to \infty }{lim}{\overline{n}}_{\ell }=\frac{{\gamma }^h}{{\gamma }^h+{\gamma }^c},$$

(A.18)

which coincides with the averaged map obtained in [1]. Note that, $\overline{n}$ together with $d,p,q\; =\; 0$ contains all the information about the final state (for the average) in the limit $\ell \; \to \; \infty$. Moreover, we have $⟨a_k^{†}{a}_k⟩\; =\; ⟨a_{-k}^{†}{a}_{-k}⟩$ (since $d\; =\; 0$), which implies in particular that the final state has a $k\; \leftrightarrow \; -k$ symmetry and, clearly, that there is a unique stationary state for the average map. From the averaged rates we can also obtain the steady-state energy

$$E_{\infty }=\; -𝜖(1\; -\; 2\overline{n})\; =\; 𝜖\frac{-{\gamma }^c+{\gamma }^h}{{\gamma }^c+{\gamma }^h}+\mathcal{\text{𝒪}}(g^2{t}^2).$$

(A.19)

This result coincides with the one found in [1], but now we show its validity for any chosen distribution of times and frequencies.

A.1.3 Standard deviation

We can also compute the standard deviation of the distribution of $n_{\ell }$ by squaring eq. (A.16), obtaining

$$n_{\ell +1}^2=n_{\ell }^2+\; 2{\eta }^2\left[|y{|}^2-|z{|}^2{n}_{\ell }\right]n_{\ell }+O({\eta }^4).$$

(A.20)

Averaging over realizations and taking the limit $\ell \; \to \; \infty$ we obtain the following variance:

$$\underset{\ell \to \infty }{lim}\overline{n_{\ell }^2}-{\overline{n_{\ell }}}^2=O({\eta }^2).$$

(A.21)

Note that the same is true of the other variables in eq. (A.14). Thus, as announced, in the weak-coupling limit each realization, i.e. the state obtained after applying $\ell \; \gg \; 1$ random cooling maps, will converge to a state close to the one obtained with the averaged map.

We emphasize that, independent of the weak-coupling limit, the standard deviation is upper bounded by $\overline{n}$, since $0\; \le n_{\ell }\le \; 1$.⁷ Thus, $\overline{n}\ll \; 1$ is sufficient in order to ensure that cooling will take place for most of the realizations. Additionally, since the average map converges to the steady state in constant time, each of the realizations will also provide cooling in constant time.

A.1.4 Cooling limit

We consider now the cooling limit, defined as ${𝜖}_{m}t\; \gg \; 1$. Let us first concentrate on the heating rate, ${\gamma }^h$. Since $|{\mathrm{\Delta }}_{m,M}+{𝜖}_M|t\; \ge {𝜖}_{m}t\; \gg \; 1$, we can expand eq. (A.17) and obtain

$${\gamma }^h=\frac{2g^2}{({\mathrm{\Delta }}_M+\; 𝜖)({\mathrm{\Delta }}_m+\; 𝜖)}+O\left(t^{-2}\right).$$

(A.22)

Here, we have taken into account that the arguments of $F$ are both positive. For ${\gamma }^c$, however, one is positive and the other negative. Assuming $({\mathrm{\Delta }}_M-𝜖)t,(𝜖\; -{\mathrm{\Delta }}_m)t\; \gg \; 1$ we obtain

$${\gamma }^c=\frac{2\pi g^{2}t}{{\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m}+O\left(1\right).$$

(A.23)

Note that the conditions above are implied by the cooling limit if we choose

$${\mathrm{\Delta }}_m={𝜖}_m∕2,{\mathrm{\Delta }}_M={𝜖}_M+{𝜖}_m∕2.$$

(A.24)

With all that, we obtain the following simple expression for $\overline{n}$ and therefore for the final state,

$$\overline{n}=\frac{2({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)}{\pi t({\mathrm{\Delta }}_M+\; 𝜖)({\mathrm{\Delta }}_m+\; 𝜖)}.$$

(A.25)

With these results, we can also compute the time required to cool the system. For that, we can use eq. (A.16) combined with the cooling limit derived here. In case we want to obtain a value of $n_{\ell }=\; O(\delta )$, where $\delta$ is a small parameter, we need to choose, according to (A.25), $t\; =\; O(1∕\delta )$. Furthermore, the total time required to reach that values is $Mt$, where $M$ is the number of cooling steps. In order to obtain that, we define $m_{\ell }=n_{\ell }-\overline{n}$, which fulfills the recursion equation $m_{\ell +1}=\; (1\; -{\gamma }^c-{\gamma }^h)m_{\ell }$. Thus imposing that $m_{\ell }=\; O(\delta )$ we need $M\; =\; O(\log <!--\; FUNCTION\; APPLICATION\; -->({\delta }^{-1})∕{\gamma }^c)$. Furthermore, since ${\gamma }^c=\; O(g^{2}t)$, and $g\; =\; O(1∕t)$ in order to respect the weak coupling limit, we have that the total time

$$T\; =\; O(Mt)\; =\; O(t^2\log <!--\; FUNCTION\; APPLICATION\; -->({\delta }^{-1}))\; =\; O({\delta }^{-2}\log <!--\; FUNCTION\; APPLICATION\; -->({\delta }^{-1})).$$

(A.26)

In case we are interested in reaching a fidelity $F\; =\; 1\; -O({\delta }_F)$ we need to choose $\delta \; ={\delta }_F∕N$, so that the total number of excitations is $\mathcal{\text{𝒪}}({\delta }_F)$. Thus, the total time scales as $T_F=\; O(N^2{\delta }_F^{-2}\log <!--\; FUNCTION\; APPLICATION\; -->(N∕{\delta }_F))$. This improves over the value obtained in [1].

A.1.5 Discrete distribution

We now consider a discrete distribution for $\mathrm{\Delta }$ composed of $R\; \propto \; t$ frequencies, equally distributed in the interval $[{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$ with a spacing $\propto t^{-1}$. We have ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}+\; (-)\delta$, with $\delta$ a constant. In the terminology used in the main text, the regime $R\; \gg \; 1$ is the many-frequency regime; in the limit $R\; \to \; \infty$ the discrete sum approaches the continuum case, derived above. We proceed as follows. Firstly, we average with respect to time. This average can then be approximated by a Lorentzian function, yielding

$${\gamma }_{\mathrm{\Delta }}^c=\frac{2g^2}{{(\mathrm{\Delta }-𝜖)}^2+{\gamma }_0^2},{\gamma }_{\mathrm{\Delta }}^h=\frac{2g^2}{{(\mathrm{\Delta }+\; 𝜖)}^2+{\gamma }_0^2},$$

(A.27)

where $2{\gamma }_0=\sqrt{6}∕t$ is the effective linewidth. In the expression for ${\gamma }_{\mathrm{\Delta }}^h$, this linewidth term can be neglected as we do not consider negative frequencies, thus $(\mathrm{\Delta }+\; 𝜖)t\; \gg \; 1$ due to the cooling limit and $\mathrm{\Delta }+\; 𝜖\; \gg {\gamma }_0$. Choosing the spacing $\mathrm{\text{x}}{\gamma }_0$, we have

$$R\; =\frac{{\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m}{\mathrm{\text{x}}{\gamma }_0},$$

(A.28)

and we find bounds for our averaged cooling rate at any given $k$:

$$min\{\underset{k_1}{\overset{c}{\gamma }},{\gamma }_{k_2}^c\}\le {\gamma }_k^c\le \; max\{\underset{k_1}{\overset{c}{\gamma }},{\gamma }_{k_2}^c\}\; =\mathcal{\text{𝒪}}(g^{2}t),$$

(A.31)

where $k_1$ is a mode located exactly in resonance with some ${\mathrm{\Delta }}_i$ and $k_2$ is located exactly between two resonance peaks. The difference ${\gamma }_{k_1}^c-{\gamma }_{k_2}^c$ tends to zero for $\mathrm{\text{x}}∕\pi \; \ll \; 1$: since $\coth <!--\; FUNCTION\; APPLICATION\; -->(y)\; -\tanh <!--\; FUNCTION\; APPLICATION\; -->(y)\; =\; 2∕\sinh <!--\; FUNCTION\; APPLICATION\; -->(2y)$, the two rates become equal for $\mathrm{\text{x}}\to \; 0$. Therefore, choosing $\mathrm{\text{x}}∕\pi \; \ll \; 1$ makes the bounds in (A.31) equal and ${\gamma }_k^c$ independent of $k$ (see figures 7 and 8). We choose $\mathrm{\text{x}}=\; 1$ as it fulfills this approximation while at the same time keeping the number of frequencies $R$ small, yielding

$${\gamma }^c\approx \frac{2\sqrt{6}\pi g^{2}t}{3({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)}.$$

(A.32)

Note that one should obtain the same result as in eq. (A.23), found in the continuum limit, but the prefactor is slightly different due to the Lorentzian approximation done in eq. (A.27). For the averaged heating rate [eq. (A.22)], the result does depend on $k$, but can be upper-bounded by ${\gamma }_{-N∕2}^h$ or roughly approximated by ${\gamma }_{-N∕4}^h$.

Let us make a few remarks. Firstly, we have chosen ${\mathrm{\Delta }}_{M,m}$ such that the modes at the edges of the spectrum can still be treated as modes in the middle of the spectrum. Choosing ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}+\; (-)\delta$, as we do for the discrete distribution, we require $\delta t\; \gg \; 1$ to remain in the cooling limit and for the analytical results to hold. However, it can be favorable to choose less frequencies, even if the edges of the spectrum are cooled unevenly and $\delta t\; \gg \; 1$ is not fulfilled. This is the case, for example, in the presence of noise, where unneeded frequencies are detrimental.

We also remark that cooling will become less effective near the phase transition. This is due to two crucial details. The first is that the choice of ${\mathrm{\Delta }}_{M,m}$ needs to be such that the edges of the spectrum ${𝜖}_{M,m}$ are surrounded by resonant frequencies on both sides. Near the phase transition, where ${𝜖}_m\to \; 0$, this cannot be achieved, since adding frequencies below ${𝜖}_m$ (possibly negative) would imply reheating other modes. The cooling limit ${𝜖}_{m}t\; \gg \; 1$ also indicates that cycle time must be inversely proportional to the gap, which will ask for an unfeasible $t$.

Finally, as mentioned before, a similar flat-rate profile was obtained in ref. [11]. However, in contrast to the choice of random cycle times and multiple frequencies proposed here, the authors of ref. [11] achieve a similar profile with time-dependent couplings.

(a)

(b)

(c)

Figure 7: Rescaled averaged cooling rate $R{\gamma }_k^c$ versus mode energy ${𝜖}_k$ for three configurations of bath frequencies. The interval $[{𝜖}_m,{𝜖}_M]$ marks the system’s spectral range; $[{\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M]$ marks the bath frequency range. (a) Too few frequencies: pronounced peaks and valleys lead to strongly non-uniform cooling. (b) Sufficiently many frequencies but short cycle times (broad linewidth ${\gamma }_0$): resonance peaks extend well beyond the spectrum. (c) Many frequencies with long cycle times (narrow linewidth): the averaged rate is approximately flat across $[{𝜖}_m,{𝜖}_M]$ and decays sharply outside, yielding the desired uniform cooling profile.

A.2 Steady-state energy in the low-noise limit

We now include depolarizing noise and derive the steady-state relative energy of mode $k$ in the low-noise limit $\kappa t\; \ll \; 1$. In the main text the noise is written as a channel [eq. (3.12)], equivalently generated by the depolarizing master equation [eq. (3.11)]. Acting on a mode pair during one cooling subcycle of duration $\tau$ it reads ${\mathcal{\text{𝒟}}}_{\tau }(\rho )\; =e^{-2\kappa \tau }\rho \; +\; (1\; -e^{-2\kappa \tau })\frac{𝟙}{4}$ and drives the average occupation towards that of the maximally mixed state, $\overline{n}=\; 1∕2$, as $\overline{n}(\tau )\; =e^{-2\kappa \tau }\overline{n}(0)\; +\; (1\; -e^{-2\kappa \tau })∕2$, i.e. $\dot{\overline{n}}=\; -2\kappa (\overline{n}-\frac{1}{2})$. In the low-noise limit $\kappa \tau \; \ll \; 1$ we expand to first order in $\kappa \tau$ and, since the depolarizing generator commutes with both the coherent evolution and the cooling Lindbladian, add this contribution directly to the single-realization recursion [eq. (A.15)],

$$n_{\ell +1}=n_{\ell }+{\eta }^2\left(|y{|}^2-|z{|}^2{n}_{\ell }\right)-\; 2\kappa \tau \left(n_{\ell }-\frac{1}{2}\right)$$

(A.33)

with corrections of order $O[{\eta }^4,∕\kappa t{)}^2,{\eta }^2\kappa ]$. Averaging with respect to $\tau$ and $\mathrm{\Delta }$ (recall $⟨\tau ⟩\; =\; t$) we obtain

This is precisely the noiseless recursion [eq. (A.16)] with the replacement ${\gamma }^{c,h}\to {\gamma }^{c,h}+\; \kappa t$, which establishes

$${\gamma }_{noisy}^{c,h}={\gamma }^{c,h}+\; \kappa t.$$

(A.35)

We can now combine these noisy rates with the results of the previous subsection to obtain the steady-state relative energy of mode $k$. From the relative energies of the individual modes one obtains the total relative energy, $e\; =\; ({\sum }_k{𝜖}_k{e}_k)∕{\sum }_k{𝜖}_k$, as a function of the cycle time $t$ and the noise rate $\kappa$; the optimization of $e$ over $t$ for any given $\kappa$ is carried out in appendix A.3.

The relative energy of mode $k$ is defined as

$$e_k=\left|\frac{E_k-E_{k,\mathrm{\text{GS}}}}{E_{k,\mathrm{\text{GS}}}}\right|,$$

(A.36)

where $E_{k,\mathrm{\text{GS}}}=\; -{𝜖}_k$ is the ground-state energy of the mode. Inserting the noiseless steady-state energy of the averaged map [eq. (A.19)] gives

$$e_k=\; 2{\overline{n}}_k+\mathcal{\text{𝒪}}(g^2{t}^2),{\overline{n}}_k=\frac{{\gamma }_k^h}{{\gamma }_k^h+{\gamma }_k^c}.$$

(A.37)

Replacing ${\gamma }_k^{c,h}\to {\gamma }_k^{c,h}+\; \kappa t$ according to eq. (A.35) gives

$$e_k=\frac{2{\gamma }_k^h+\; 2\kappa t}{{\gamma }_k^h+{\gamma }_k^c+\; 2\kappa t},$$

(A.38)

which is the expression quoted in the main text [cf. eq. (4.10)]. Now, we can choose any distribution of frequencies and calculate the rates.

A.2.1 Discrete distribution

We now insert the explicit averaged heating and cooling rates derived in the weak-coupling and cooling limits [eqs. (A.22) and (A.32)] for the discrete distribution (see appendix A.1.5). Dividing the numerator and denominator of eq. (A.38) by $2g^2$, we obtain

To make the structure of eq. (A.39) transparent, we collect the three contributions into the coefficients

where $A_k$ originates from the heating rate [eq. (A.22)], $Ct$ from the cooling rate [eq. (A.32)], and $Bt$ from the noise term [eq. (A.35)]. The relative energy in eq. (A.39) then can be written as

From eq. (A.41) one clearly sees that cooling ($e_k\to \; 0$) is obtained when $Ct\; \gg A_k,Bt$, which translates into $t\; \gg \; ({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)∕({\mathrm{\Delta }}_m{\mathrm{\Delta }}_M)$ and $\kappa ∕g^2\ll \; 1∕({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)$. Taking ${\mathrm{\Delta }}_m,{\mathrm{\Delta }}_M$ as in the choice ${\mathrm{\Delta }}_m={𝜖}_m∕2,{\mathrm{\Delta }}_M={𝜖}_M+{𝜖}_m∕2$ used above, the first condition is nothing but the cooling limit ${𝜖}_{m}t\; \gg \; 1$, whereas the second, combined with $\kappa \; \ll {𝜖}_m$, implies $\kappa \; \ll \; g$. Thus we can only expect cooling under this condition. In addition, if we merely impose $Bt\; \gg A_k$, i.e. $\kappa t∕g^2\gg \; 1∕({\mathrm{\Delta }}_m{\mathrm{\Delta }}_M)$, we obtain a simple expression for the final energy:

$$e^{cool}\approx \frac{1}{1\; +\frac{\sqrt{2}\pi g^2}{\sqrt{3}\kappa ({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)}}.$$

(A.42)

(a)

(b)

(c)

Figure 9: Relative energy $e$ versus the ratio $\kappa ∕g^2$ for cooling with three cycle times: (a) $t\; =\; 10$, (b) $t\; =\; 50$, (c) $t\; =\; 100$, and four coupling strengths $g\; =\; 0.003$ (blue), $0.005$ (orange), $0.013$ (green), $0.025$ (purple). The dashed black line shows the analytic prediction from eq. (A.42). Agreement is good when $\kappa ∕g^2\gtrsim \; 1$ and $gt\; \lesssim \; 0.5$; deviations at smaller $\kappa ∕g^2$ or larger $gt$ indicate breakdown of the long-time or weak-coupling approximations. Parameters: $N\; =\; 20$, $𝜃\; =\; \pi ∕8$, ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}\pm \; 0.1$.

A.3 Optimal cycle time and energy

Having identified the regime in which cooling occurs, we now optimize the average cycle time $t$ at fixed noise rate $\kappa$ to obtain the lowest achievable relative energy. The parameter regime in which eq. (A.42) approximates the relative energy well is illustrated in figure 9: we can see that the weak-coupling condition prevents us from choosing $t$ arbitrarily large. Thus, we will now consider a constrained optimization of $e_k$ in which the product $gt$ is kept fixed and denote, as before,

Under this constraint $g\; =\; \eta ∕t$, so the noise term in eq. (A.41) grows as $Bt\; =\; \kappa {\eta }^{-2}{t}^3$. Substituting into eq. (A.39) gives

Unlike the unconstrained case, the noise now grows as $t^3$ while the cooling grows linearly in $t$, so $e_k(t)$ has a minimum at finite $t$. Setting $\partial e_k∕\partial t\; =\; 0$ yields

This can be solved exactly, but the scaling is clearer under a mild approximation: since $A_k,C\; =\mathcal{\text{𝒪}}(1)$ and $Ct\; \gg A_k$ in the cooling limit, the cubic term dominates the left-hand side, giving $t^3\approx A_k{\eta }^2∕\kappa$. The optimal cycle time and energy of mode $k$ are then

These expressions hold provided $\kappa t_{opt,k}^{cool}\ll \; 1$, consistent with the low-noise assumption. This scaling is expected to hold for the whole system, even though the actual optimal times differ slightly from mode to mode (since $A_k$ depends on $k$). We can also observe that, since the cooling rate increases with the value of $gt$, it is advisable to choose the highest possible value of $\eta$ that still satisfies the weak coupling condition in order to achieve the best efficiency. For all numerical results and figures, we use $\eta \; =\; 0.5$, as it numerically proves to be close to the limit where the weak-coupling assumption does not hold anymore, but still holds. We further approximate $A_k\approx A_{N∕4}$, which acts as a representative value across all $A_k$, to obtain an approximation for the total relative energy $e_{opt}^{cool}$. This approximation is numerically verified to be close to the exact value of $e_{opt}^{cool}$ for the whole spectrum and range of $\kappa$ studied (see figure 4).

As a final note, in the analysis above we have treated $R\; \gg \; 1$ as a continuous variable to make the functions well-behaved, but in actual simulations $R\; \in \; \mathbb{N}$, which produces discrete jumps whenever a new frequency is added. In special cases like vanishing spectrum or high noise rates, using only one or a few frequencies might be advantageous over the $R\; \gg \; 1$ regime. We discuss this in the next subsection.

A.4 The few frequency case

In the previous subsections, we have analyzed the many-frequency regime, $R\; \gg \; 1$, which allowed us to cool all modes in the system evenly, while at the same time keeping the effect of noise low. We now consider two alternatives: the single-frequency ($R\; =\; 1$) and few-frequency ($R\; \gtrsim \; 1$) regimes. These regimes are advantageous when the spectrum is narrow enough that one or only a few resonant frequencies cool all modes, or when the noise rate is large enough that the optimized cycle time is short. Adding unnecessary frequencies should be avoided, because it increases the runtime of the algorithm and makes the impact of noise larger. We will analyze in this section where the $R\; =\; 1$ regime can prove advantageous, and what corrections can be made to our algorithm in order to bridge the gap between the $R\; \gg \; 1$ and the $R\; =\; 1$ regimes. Finally, we find a formula for the scaling with noise rate $\kappa$ that can be applied in the intermediate regime.

We will begin by analyzing an important region of our model: the edge of the trivial phase, $𝜃\; =\; \pi ∕2$. In this case, only a single frequency is required, since all $k$-mode energies ${𝜖}_k=\; 1$ coincide. Therefore, we have $R\; =\; 1$ independent of the cycle time $t$ and of the noise rate $\kappa$. This causes the dependence of the steady state energy on $\kappa$ to differ from that observed in the many-frequency regime. Following similar steps as in the previous subsections, we obtain

This protocol yields a steady state ${\rho }_k={\rho }_{00}|0⟩⟨0|+\; (1\; -{\rho }_{00})|1⟩⟨1|$,⁸ that is the same for every mode, with

$${\rho }_{00}\approx \; 1\; -\frac{D}{2}{\kappa }^{2∕3}.$$

(A.51)

The scaling in the single-frequency scenario is $e\; \sim {\kappa }^{2∕3}$, clearly different from the many-frequency case detailed previously. Now, we will show a way to bridge the area between the two regimes. To that end, we introduce the following corrections to the protocol:

• $R$ is taken to be a natural number: $R\; =\; \lfloor \sqrt{2∕3}({\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m)t\rfloor \; +\; 1$, so that when ${\mathrm{\Delta }}_M-{\mathrm{\Delta }}_m\to \; 0$ the single-frequency result is recovered, and in the limit $R\; \gg \; 1$ one extra frequency has a negligible effect.

• We take ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}+\; (-)f(𝜃)$, where $f(𝜃)$ is a positive function that is roughly constant for most values of $𝜃$ (ensuring that the full spectrum is still covered as explained in previous subsections) but goes to zero in the limits $𝜃\; \to \; 0$ and $𝜃\; \to \pi ∕2$. We can choose, for example, the following approximation to a square pulse:

  $$f(𝜃)\; =\frac{\tanh <!--\; FUNCTION\; APPLICATION\; -->(20𝜃\; -\; 1)\; -\tanh <!--\; FUNCTION\; APPLICATION\; -->(20𝜃\; -\; 10\pi \; +\; 1)}{20}.$$

  (A.52)

  In general, for this choice of $f(𝜃)$, ${\mathrm{\Delta }}_{M(m)}\approx {𝜖}_{M(m)}\pm \; 0.1$ for most values of $𝜃$. But for $𝜃\; \in \; [0,0.1]$ (and analogously for $𝜃\; \in \; [\pi ∕2\; -\; 0.1,\pi ∕2]$) the function behaves as a linear decay: $f(𝜃)\; \approx 𝜃$.

In the regime $𝜃\; \in \; [0,0.1]$, $R$ contains both a term linear in $t$ and a constant contribution (the additional $+\; 1$), neither of which can be neglected. Furthermore, the dependence on $t$ is not smooth, as $R$ exhibits discrete jumps due to its integer nature; consequently, there exists a range of cycle times that yield the same value of $R$. In addition, the coefficients ${\gamma }_k^{c,h}$ depend on the number of frequencies and therefore also display discontinuous behavior. Taking these features into account, the determination of the optimal cycle time and steady-state energy in this regime is most reliably carried out numerically.

Moreover, since the optimal value of $t$ decreases as $\kappa$ increases, a crossover from ${\kappa }^{1∕3}$ to ${\kappa }^{2∕3}$ scaling can occur for sufficiently large $\kappa$ across the spectrum, as illustrated in figure 10. Finding the true analytical expression for the few-frequency regime between the $\mathcal{\text{𝒪}}({\kappa }^{1∕3})$ and $\mathcal{\text{𝒪}}({\kappa }^{2∕3})$ regimes would require summing over all different $k$ modes, each with their own $k$-dependent cooling and heating rates (see figure 7b). However, the total relative energy can be well approximated by simply combining the contributions from both analytically tractable limits:

$$e_{opt}^{cool}\approx D{\kappa }^{2∕3}+\frac{3}{C}{(A^2{\eta }^{-2}\kappa )}^{1∕3}.$$

(A.53)

This result is supported by the exact numerics shown in figure 10.

Note that the number of additional frequencies could, in principle, also be optimized, suggesting the existence of an optimal functional form of $f(𝜃)$; however, identifying this form lies beyond the scope of the present work.

(a)

(b)

(c)

(d)

Figure 10: Optimal relative energy $e_{opt}^{cool}$ versus noise strength $\kappa$ (both log scale) for $N\; =\; 100$, $gt\; =\; 0.5$, and ${\mathrm{\Delta }}_{M(m)}={𝜖}_{M(m)}\pm f(𝜃)$ in our cooling protocol. Blue circles: exact numerical optimization over $t$. Dashed black lines: combined approximation from eq. (A.53). Dashed green lines: multifrequency scaling $\propto {\kappa }^{1∕3}$. Dashed orange lines: single-frequency scaling $\propto {\kappa }^{2∕3}$. (a) $𝜃\; =\; \pi ∕1000$: the single-frequency term dominates. (b) $𝜃\; =\; \pi ∕100$: both terms contribute. (c) $𝜃\; =\; \pi ∕50$: multifrequency term begins to dominate. (d) $𝜃\; =\; \pi ∕20$: the ${\kappa }^{1∕3}$ many-frequency scaling holds over most of the range.

B Analytic derivation of the adiabatic evolution

In this appendix we derive how the relative energy of the adiabatic algorithm scales with the noise strength $\kappa$. We use here the well-known single-mode adiabatic perturbation theory and Landau-Zener formulas [3, 40, 41] and apply them to our specific model. We then sum over all momentum modes to obtain the total relative energy, and optimize over the evolution time $T$ in the presence of depolarizing noise. We find formulas for the relative energy in the trivial phase (appendix B.1) and the topological phase (appendix B.2), then optimize the relative energy in both phases taking into account the effect of noise (appendix B.3), and finish by modifying our formulas in the case of a realistic initial state (appendix B.4).

B.1 Trivial phase: first-order adiabatic perturbation theory

Let us first review the adiabatic approximation and then apply it to our protocol. We start with the Schrödinger equation for a time-dependent two-level system corresponding to a single momentum mode $k$. For clarity we drop the $k$ subscript and write

$$i\frac{d}{dt}|\psi (t)⟩=\; H(t)|\psi (t)⟩,$$

(B.1)

with

$$H(t)\; =\left(\begin{array}{cc}\hfill f(t)\hfill & \hfill r(t)\hfill \\ \hfill r(t)\hfill & \hfill -f(t)\hfill \end{array}\right),$$

(B.2)

and $t\; \in \; [0,T]$. We use the rescaled variable $s\; =\; t∕T\; \in \; [0,1]$. The eigenvalues and eigenvectors (in terms of $s$) are:

We expand $|\psi (t)⟩=c_0(t)|0(t)⟩+c_1(t)|1(t)⟩$ in the instantaneous eigenbasis. Substituting into the Schrödinger equation and projecting onto the instantaneous eigenbasis gives:

The adiabatic approximation consists of neglecting the off-diagonal coupling $c_1⟨0|\dot{1}⟩$ in the first equation, which is valid when the spectral gap is large. Solving for $c_0(0)\; =\; 1$ we obtain

Substituting back into the second equation and solving for $c_1(0)\; =\; 0$ yields

Using $s\; =\; t∕T$ and denoting $v\; =\; 1∕T$ we have

$$c_1(s)\; ={\int }_0^s{e}^{-i{\gamma }_{10}(s^{\prime })}{e}^{\frac{i}{v}{\omega }_{10}(s^{\prime })}⟨1(s^{\prime })|\dot{0}(s^{\prime })⟩\mathrm{\text{d}}{s}^{\prime },$$

(B.11)

with Berry and dynamical phases:

where $n,m$ take the values $0,1$. To extract the leading behavior for large $T$ (small $v$) in eq. (B.11), we use the following identity for integrals of the form

$$I\; ={\int }_0^s\mathrm{\text{d}}{s}^{\prime }B(s^{\prime })e^{\frac{1}{v}{\int }_0^{s^{\prime }}\mathrm{\text{d}}{s}^{″}C(s^{″})},$$

(B.14)

where we can expand them in growing terms of $v$ by iteratively applying the identity

The first term gives us the first order in $v$; the second gives us an integral of the same type as $I$ but with a factor of $v$ already in front and $B$ replaced by $\frac{d}{\mathrm{\text{d}}{s}^{\prime }}\left(\frac{B(s^{\prime })}{C(s^{\prime })}\right)$, therefore contributing only to higher orders of $v$. We keep only the first order in $v$.

Let us now apply this to our choice of $H(s)$ (see eq. (4.1)). The Berry phase vanishes: ${\gamma }_{10}(s)\; =\; 0$, so eq. (B.11) has the shape of eq. (B.14). The resulting transition probability to the excited state after the complete evolution is, up to first order in $v$,

$$|c_1(1){|}^2\approx \frac{1}{T^2}{\left|e^{iT{\omega }_{10}(1)}\frac{⟨1(1)|\dot{0}(1)⟩}{2E_1(1)}-\frac{⟨1(0)|\dot{0}(0)⟩}{2E_1(0)}\right|}^2.$$

(B.16)

We now map this to our specific Hamiltonian. Starting from ${𝜃}_i=\; \pi ∕2$ and targeting ${𝜃}_f$ via the linear interpolation $H(s)\; =H_0(1\; -s)\; +H_{f}s$, the matrix elements are

where ${\varphi }_k=\; 2\pi k∕N$.

The dynamical phase integral evaluates to

with

The matrix elements in eq. (B.16) for $s\; >\; 0$ are given by:

$$⟨1(s)|\dot{0}(s)⟩=\frac{r}{s}\left(1\; +\frac{f}{\sqrt{f^2+r^2}}\right)\frac{1}{r^2+{(f\; +\sqrt{f^2+r^2})}^2},$$

(B.22)

whereas for $s\; =\; 0$ we have

$$⟨1(0)|\dot{0}(0)⟩=\frac{1}{2}\cos <!--\; FUNCTION\; APPLICATION\; -->{𝜃}_f\sin <!--\; FUNCTION\; APPLICATION\; -->{\varphi }_k.$$

(B.23)

Since we are interested in scalings rather than exact prefactors, we do not evaluate these expressions further. Summing over all modes, the total energy after the adiabatic evolution in the absence of noise is

$$E\; =\underset{k}{\sum }-{𝜖}_k\left(1\; -|c_1(k){|}^2\right)=E_{GS}\left(1\; -\frac{Ã}{T^2}\right).$$

(B.24)

The prefactor $Ã$ can be calculated numerically. We have assumed that the time dependence inside the absolute value in eq. (B.16) can be neglected, since it is only a phase and does not significantly alter the value of $|c_1(k){|}^2$; consequently, $Ã$ is approximately time-independent. One can further approximate $Ã\; \approx T^2|c_1(N∕4){|}^2\approx \left({\cos <!--\; FUNCTION\; APPLICATION\; -->}^2{𝜃}_f\right)∕8$, which can be verified numerically.

B.2 Topological phase: Landau-Zener transition

When the adiabatic path crosses the phase transition, the perturbative treatment fails for modes in the vicinity of the gap-closing point. In this regime, the energy can instead be determined using the Landau-Zener formalism [3, 40, 41]. We apply the standard Landau-Zener formulas to our adiabatic evolution.

This approximation is most accurate when the energy gaps in the initial and final configurations are large. To extract an analytic approximation, we follow a path starting from $𝜃\; =\; \pi ∕2$ and going to $𝜃\; =\; 0$. We redefine our linear path such that $H_k(s)\; =H_k^{\pi ∕2}(0.5\; -s)\; +H_k^0(0.5\; +\; s)$, where $s\; =\; t∕T$ and we run $t$ from $-\; T∕2$ to $T∕2$. Close to the transition point $s\; \approx \; 0$ and close to the mode where the gap closes, ${\delta }_k={\varphi }_k-\pi \; \approx \; 0$ (with $k\; \in \; [0,N]$), the eigenenergies of the instantaneous Hamiltonian can be approximated by

$${𝜖}_k({\delta }_k,t)\; \approx \sqrt{\frac{4t^2}{T^2}+\frac{{\delta }_k^2}{4}}.$$

(B.25)

This dispersion relation can be identified with that of the standard Landau-Zener Hamiltonian

which, by comparison, gives

From the standard Landau-Zener result, the transition probability is (with $\hslash \; =\; 1$):

$$P_{LZ}^k=\exp <!--\; FUNCTION\; APPLICATION\; -->\left(-\frac{\pi {\mathrm{\Delta }}_{LZ}^2}{v_{LZ}}\right)=\exp <!--\; FUNCTION\; APPLICATION\; -->\left(-\frac{\pi }{8}T{\delta }_k^2\right).$$

(B.30)

This gives a relative energy per mode $e_k=\; 2P_{LZ}^k$. As expected, the smaller the gap of the mode, the longer the evolution time required to achieve a low excitation probability. A true closing of the gap occurs at $k\; =\; N∕2$ (i.e., ${\delta }_k=\; 0$), but in the limit $N\; \gg \; 1$ this single-mode contribution becomes negligible.

This approximation does not depend on the target angle ${𝜃}_f$ and was derived for the specific path to ${𝜃}_f=\; 0$ as an illustration.

A result for general ${𝜃}_f$ can be extracted in a similar fashion:

We recover $G({𝜃}_f=\; 0)\; =\; 1∕8$, but for other values of ${𝜃}_f$ in the topological phase, $G({𝜃}_f)$ changes significantly, e.g., $G(\pi ∕8)\; \approx \; 0.23$.

The derivation of the total energy is more involved and, in general, does not yield a closed-form expression. Following the Kibble-Zurek argument [22, 23, 36–38], for a large enough system size the energy can be approximated by

where the integral has been extended to infinity since the contribution of modes with large gaps is negligible. The value of the integral depends on the target angle ${𝜃}_f$. For ${𝜃}_f=\; 0$ or ${𝜃}_f=\; \pi ∕2$, one has ${𝜖}_k=\; 1$ for all $k$ and a closed form exists; in the general case, further approximation is needed. For large $T$, the Gaussian factor is sharply peaked in $x$ and the energy factor is approximately constant over the range where the Gaussian is concentrated. Evaluating the energy at the peak of the Gaussian and taking it out of the integral gives

To obtain a closed equation for the relative energy $e$, an approximation for $E_{GS}$ is required, which we derive below.

Approximating $E_{GS}$. In the large system size limit, the ground state energy is given by the following integral:

$$E_{GS}(𝜃)\; =\; -\frac{N}{\pi }{\int }_0^{\pi }\sqrt{1\; +\sin <!--\; FUNCTION\; APPLICATION\; -->(2𝜃)\cos <!--\; FUNCTION\; APPLICATION\; -->x}\mathrm{\text{d}}x.$$

(B.35)

To approximate this integral, we construct a function that has the same expansion as the $|\sin <!--\; FUNCTION\; APPLICATION\; -->(2𝜃)|\; \ll \; 1$ approximation of the integral, is symmetric with respect to $𝜃\; =\; \pi ∕4$, and lies close to the exact value for $𝜃\; =\; \pi ∕4$:

$$E_{GS}(𝜃)\; \approx -N\sqrt{1\; -\frac{{\sin <!--\; FUNCTION\; APPLICATION\; -->}^2(2𝜃)}{8}-\frac{{\sin <!--\; FUNCTION\; APPLICATION\; -->}^4(2𝜃)}{17}}.$$

(B.36)

This approximation can be numerically verified to be within $1\%$ of the true value of the integral, and allows us to approximate the ground-state energy for any given $𝜃$ in the large $N$ limit. It can also be verified to lie close to the Padé-Chebyshev approximant for the same function.

B.3 Adding noise

(a)

(b)

Figure 11: Relative energy $e$ versus total evolution time $T$ (log-log scale) for adiabatic evolution targeting (a) ${𝜃}_f=\; \pi ∕8$ and (b) ${𝜃}_f=\; \pi ∕3$ at several noise rates $\kappa$, for system size $N\; =\; 20$. Solid lines: exact numerical simulation. Dashed lines: analytic predictions. In the presence of noise, each curve exhibits a minimum at an optimal $T$, reflecting the trade-off between non-adiabatic excitations and depolarizing noise. For larger $T$, the energy increases as noise effects dominate. The analytical predictions accurately capture both the scaling behavior and the position of these minima.

For depolarizing noise of strength $\kappa$, the noise channel commutes with the unitary evolution (see ref. [1] and section 3.3). For our model, the maximally mixed state has $E\; =\; 0$. Then, the energy in the presence of noise is related to the energy in the noiseless setting by an overall damping factor as

$$E_{\text{noisy}}=e^{-2\kappa T}{E}_{\text{ideal}}.$$

(B.37)

Figure 11 illustrates the effect of noise on the adiabatic evolution targeting ${𝜃}_f=\; \pi ∕3$ (within the trivial phase, no gap closing). It shows that the analytical predictions accurately capture both the scaling behavior and the position of these minima. For short $T$, the energy decreases as $T^{-2}$ and is dominated by non-adiabatic excitations present even in the absence of noise. For large $T$, depolarization effects become dominant. The optimal $T$ arises from balancing these two contributions.

We can now determine the scalings for the adiabatic algorithm in the presence of noise, for both cases (crossing or not crossing the phase transition).

Gapped regime (no gap closing, ${𝜃}_f>\; \pi ∕4$). When the path does not cross the phase transition and for sufficiently small noise rate $\kappa$, the total energy can be easily determined by multiplying eq. (B.24) with the exponential damping factor resulting from the depolarizing noise. We obtain

$$E\; =e^{-2\kappa T}\left(1\; -\frac{Ã}{T^2}\right)E_{GS}\approx E_{GS}\left(1\; -\; 2\kappa T\; -\frac{Ã}{T^2}\right),$$

(B.38)

where the approximation results from an expansion around small $\kappa T$. Figure 11 illustrates the agreement between analytics and numerics in this scenario. Differentiating with respect to $T$ and setting the derivative to zero yields the optimal time and energy:

where $e_{opt}$ is calculated using eq. (2.1).

This result is consistent with the expectation that, in the absence of noise, the adiabatic algorithm in the limit $T\; \to \; \infty$ produces the exact ground state.

Gapless regime (gap closing, ). The analysis when the path crosses the phase transition requires more care, since each mode contributes differently to the total energy, yet all modes evolve for the same time $T$. We assume large $N$ and large $T$, multiply eq. (B.34) with the exponential damping factor resulting from the depolarizing noise, and use the same approximation as before to obtain:

with

In order to ease the determination of the optimal $T$ we approximate this expression further. Using that $\tilde{B}=\mathcal{\text{𝒪}}(1)$ and $\kappa T\; \ll \; 1$, we have that the term $2\kappa \tilde{B}\sqrt{T}$ is negligible compared to $\tilde{B}∕\sqrt{T}$ and we can simplify eq. (B.42):

Differentiating with respect to $T$ and optimizing yields:

These scalings are consistent with the numerical results presented in the main text (figure 4). The prefactor $\tilde{B}$ is monotonically decreasing from $𝜃\; =\; 0$, where $\tilde{B}\approx \; 0.9$, to $𝜃\; =\; \pi ∕4$, where $\tilde{B}=\; 0$. This is consistent with the intuition that the deeper into the topological phase the target lies, the more difficult it becomes to prepare the ground state via the adiabatic algorithm, as the effects of finite evolution time cause some modes to have non-zero transition probability to their excited states.

Effects of finite system size. In figure 4, for adiabatic trajectories that cross the critical point, the relative energy exhibits a crossover from ${\kappa }^{1∕3}$ scaling to a constant value independent of $\kappa$ for sufficiently small $\kappa$. This constant contribution arises from the mode $k\; =\; N∕2$, whose gap closes at the phase transition. The excitation of this mode is unavoidable when crossing the critical point, irrespective of $T$.

As a result, in the limit $\kappa \; =\; 0$, the optimal evolution time $T\; \to \; \infty$yields an energy $E\; =E_{GS}+\; 2{𝜖}_{N∕2}$. The crossover between the two regimes occurs at ${\kappa }^{1∕3}\sim {𝜖}_{N∕2}∕|E_{GS}|$, as confirmed numerically in figure 5. For larger system sizes, this crossover shifts to smaller values of $\kappa$, and vanishes in the thermodynamic limit $N\; \gg \; 1$, where the contribution of a single mode becomes negligible.

B.4 Realistic adiabatic evolution

Instead of starting from the ideal product state ${|0⟩}^{\otimes N}$, we consider a more realistic scenario in which the initial state is already affected by noise. To this end, we apply the cooling procedure at $𝜃\; =\; \pi ∕2$ for a given noise rate $\kappa$ (see appendix A.4). This yields a state ${\rho }^{cool}(\pi ∕2)$ that is diagonal in the mode basis, with each mode in a mixture ${\rho }_k={\rho }_{00}|0⟩⟨0|+\; (1\; -{\rho }_{00})|1⟩⟨1|$ and ${\rho }_{00}\approx \; 1\; -\; (D∕2){\kappa }^{2∕3}$.

Since all derivations in the previous subsections for an initial state $|0⟩$ apply analogously to an initial state $|1⟩$ (due to the symmetry of the spectrum), we can apply the derivations independently to both cases, weight the resulting energies by ${\rho }_{00}$ and $1\; -{\rho }_{00}$, respectively, and sum them to obtain the total energy:

where the factor $D$ has been defined in eq. (A.50) (see appendix A.4 for the full derivation). Crucially, this does not alter the optimal times, so the optimizations performed for the ideal initial state remain valid. The final energies change in different ways depending on the target phase. In the trivial phase, which already has a ${\kappa }^{2∕3}$ scaling, the correction merely modifies the prefactor. In the topological phase, the imperfect initialization contributes only a subleading additive term leading to

In the topological phase, the correction from the imperfect initial state ($D{\kappa }^{2∕3}$) becomes negligible for sufficiently small $\kappa$: this is consistent with the fact that, at low noise rates, cooling can prepare the initial state very close to ${|0⟩}^{\otimes N}$. For sufficiently large $\kappa$, however, both contributions are significant and the curve deviates from the ideal ${\kappa }^{1∕3}$ behavior.

References

[1]

D. Molpeceres, S. Lu, J.I. Cirac and B. Kraus, Quantum algorithms for cooling: A simple case study, Phys. Rev. Res. 7 (2025) 033162 [arXiv:2503.24330] [INSPIRE].

[2]

E. Farhi, J. Goldstone and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv:1411.4028 [INSPIRE].

[3]

T. Albash and D.A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90 (2018) 015002 [arXiv:1611.04471] [INSPIRE].

[4]

L. Zhou et al., Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Phys. Rev. X 10 (2020) 021067 [arXiv:1812.01041] [INSPIRE].

[5]

M. Cerezo et al., Variational quantum algorithms, Nature Rev. Phys. 3 (2021) 625 [arXiv:2012.09265] [INSPIRE].

[6]

B. Kraus et al., Preparation of entangled states by quantum Markov processes, Phys. Rev. A 78 (2008) 042307 [arXiv:0803.1463] [INSPIRE].

[7]

F. Verstraete, J.I. Cirac and J.I. Latorre, Quantum circuits for strongly correlated quantum systems, Phys. Rev. A 79 (2009) 032316 [arXiv:0804.1888] [INSPIRE].

[8]

S. Diehl et al., Quantum states and phases in driven open quantum systems with cold atoms, Nature Phys. 4 (2008) 878 [arXiv:0803.1482] [INSPIRE].

[9]

B.M. Terhal and D.P. DiVincenzo, On the problem of equilibration and the computation of correlation functions on a quantum computer, Phys. Rev. A 61 (2000) 22301 [quant-ph/9810063] [INSPIRE].

[10]

J. Lloyd et al., Quasiparticle Cooling Algorithms for Quantum Many-Body State Preparation, PRX Quantum 6 (2025) 010361 [arXiv:2404.12175] [INSPIRE].

[11]

J. Lloyd and D.A. Abanin, Quantum thermal state preparation for near-term quantum processors, arXiv:2506.21318 [INSPIRE].

[12]

Y. Zhan et al., Rapid Quantum Ground State Preparation via Dissipative Dynamics, Phys. Rev. X 16 (2026) 011004 [arXiv:2503.15827] [INSPIRE].

[13]

Z. Ding, Y. Zhan, J. Preskill and L. Lin, Simple and efficient end-to-end quantum thermal and ground state preparation, arXiv:2508.05703 [DOI:10.1038/s41567-026-03389-y] [INSPIRE].

[14]

X. Mi et al., Stable quantum-correlated many-body states through engineered dissipation, Science 383 (2024) adh9932 [arXiv:2304.13878] [INSPIRE].

[15]

D. Hahn, S.A. Parameswaran and B. Placke, Efficient quantum thermal state preparation via local driving: Lindbladian simulation with provable guarantees, Phys. Rev. B 114 (2026) 014302 [arXiv:2505.22816] [INSPIRE].

[16]

J. Langbehn et al., Universal cooling of quantum systems via randomized measurements, arXiv:2506.11964 [INSPIRE].

[17]

C. Ramon-Escandell, A. Prositto and D. Segal, Thermal state preparation by repeated interactions at and beyond the Lindblad limit, Phys. Rev. Res. 7 (2025) 043012 [arXiv:2506.12166] [INSPIRE].

[18]

K. Wang and Z. Ding, Beyond Lindblad Dynamics: Rigorous Guarantees for Thermal and Ground State Preservation under System Bath Interactions, arXiv:2512.03457 [INSPIRE].

[19]

D. Motlagh, M.S. Zini, J.M. Arrazola and N. Wiebe, Ground State Preparation via Dynamical Cooling, arXiv:2404.05810 [INSPIRE].

[20]

C.-F. Chen, M. Kastoryano, F.G.S.L. Brandão and A. Gilyén, Efficient quantum thermal simulation, Nature 646 (2025) 561 [INSPIRE].

[21]

H. Schlömer, L. Chen, S.F. Yelin and H.-Y. Hu, Model-agnostic cooling algorithms for strongly interacting fermions, arXiv:2605.01029 [INSPIRE].

[22]

W.H. Zurek, U. Dorner and P. Zoller, Dynamics of a quantum phase transition, Phys. Rev. Lett. 95 (2005) 105701 [cond-mat/0503511] [INSPIRE].

[23]

J. Dziarmaga, Dynamics of a Quantum Phase Transition: Exact Solution of the Quantum Ising Model, Phys. Rev. Lett. 95 (2005) 245701 [INSPIRE].

[24]

J.R. McClean et al., Barren plateaus in quantum neural network training landscapes, Nature Commun. 9 (2018) 4812 [arXiv:1803.11173] [INSPIRE].

[25]

J.F. Poyatos, J.I. Cirac and P. Zoller, Quantum Reservoir Engineering with Laser Cooled Trapped Ions, Phys. Rev. Lett. 77 (1996) 4728 [INSPIRE].

[26]

Z. Ding, C.-F. Chen and L. Lin, Single-ancilla ground state preparation via Lindbladians, Phys. Rev. Res. 6 (2024) 033147 [arXiv:2308.15676] [INSPIRE].

[27]

A.M. Childs, E. Farhi and J. Preskill, Robustness of adiabatic quantum computation, Phys. Rev. A 65 (2002) 012322 [quant-ph/0108048] [INSPIRE].

[28]

M. Sarovar and K.C. Young, Error suppression and error correction in adiabatic quantum computation II: non-equilibrium dynamics, New J. Phys. 15 (2013) 125032 [arXiv:1307.5892] [INSPIRE].

[29]

B.F. Schiffer, A.F. Rubio, R. Trivedi and J.I. Cirac, The quantum adiabatic algorithm suppresses the proliferation of errors, arXiv:2404.15397 [INSPIRE].

[30]

P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals Phys. 57 (1970) 79 [INSPIRE].

[31]

M. Greiter, V. Schnells and R. Thomale, The 1D Ising model and the topological phase of the Kitaev chain, Annals Phys. 351 (2014) 1026.

[32]

G.B. Mbeng, A. Russomanno and G.E. Santoro, The quantum Ising chain for beginners, SciPost Phys. Lect. Notes 82 (2024) 1 [arXiv:2009.09208] [INSPIRE].

[33]

A. Kitaev, Unpaired Majorana fermions in quantum wires, Phys. Usp. 44 (2001) 131 [cond-mat/0010440] [INSPIRE].

[34]

G. Ortiz et al., Many-Body Characterization of Particle-Conserving Topological Superfluids, Phys. Rev. Lett. 113 (2014) 267002 [arXiv:1407.3793] [INSPIRE].

[35]

V. Murg and J.I. Cirac, Adiabatic Time Evolution in Spin-Systems, Phys. Rev. A 69 (2004) 042320 [quant-ph/0309026] [INSPIRE].

[36]

T.W.B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9 (1976) 1387 [INSPIRE].

[37]

W.H. Zurek, Cosmological Experiments in Superfluid Helium?, Nature 317 (1985) 505 [INSPIRE].

[38]

A. Polkovnikov, Universal adiabatic dynamics in the vicinity of a quantum critical point, Phys. Rev. B 72 (2005) 161201 [INSPIRE].

[39]

A. del Campo and W.H. Zurek, Universality of phase transition dynamics: Topological Defects from Symmetry Breaking, Int. J. Mod. Phys. A 29 (2014) 1430018 [arXiv:1310.1600] [INSPIRE].

[40]

L.D. Landau, A theory of energy transfer. 2, Phys. Z. Sowjetunion 2 (1932) 46 [INSPIRE].

[41]

C. Zener, Nonadiabatic crossing of energy levels, Proc. Roy. Soc. Lond. A 137 (1932) 696 [INSPIRE].

[42]

J.J. Sakurai and J. Napolitano, Modern Quantum Mechanics, 3 ed., Cambridge University Press (2020) [DOI:10.1017/9781108587280] [INSPIRE].

[43]

X. Chen, Z.-C. Gu, Z.-X. Liu and X.-G. Wen, Symmetry-Protected Topological Orders in Interacting Bosonic Systems, Science 338 (2012) 1604 [arXiv:1301.0861] [INSPIRE].

[44]

S. Bravyi, M.B. Hastings and F. Verstraete, Lieb-Robinson Bounds and the Generation of Correlations and Topological Quantum Order, Phys. Rev. Lett. 97 (2006) 050401 [quant-ph/0603121] [INSPIRE].

[45]

X. Chen, Z.C. Gu and X.G. Wen, Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order, Phys. Rev. B 82 (2010) 155138 [arXiv:1004.3835] [INSPIRE].

[46]

D. Guéry-Odelin et al., Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys. 91 (2019) 045001 [arXiv:1904.08448] [INSPIRE].

[47]

A.M. Childs et al., Theory of Trotter Error with Commutator Scaling, Phys. Rev. X 11 (2021) 011020 [arXiv:1912.08854] [INSPIRE].

[48]

J.T. Barreiro et al., An open-system quantum simulator with trapped ions, Nature 470 (2011) 486 [arXiv:1104.1146] [INSPIRE].

[49]

G. Kishony, M.S. Rudner, A. Rosch and E. Berg, Gauged Cooling of Topological Excitations and Emergent Fermions on Quantum Simulators, Phys. Rev. Lett. 134 (2025) 086503 [arXiv:2310.16082] [INSPIRE].

[50]

G. Kishony, M.S. Rudner and E. Berg, Efficiently preparing chiral states via fermionic cooling on bosonic quantum hardware, Commun. Phys. 8 (2025) 73 [arXiv:2409.02177] [INSPIRE].

[51]

R. Wang et al., From exponential to quadratic: optimal control for a frustrated Ising ring model, Quantum Sci. Technol. 10 (2025) 035052 [arXiv:2502.17015] [INSPIRE].

[52]

K. Thengil, V.R. Arezzo and G.E. Santoro, Digital techniques for the frustrated Ising ring: the role of counter-diabatic terms, arXiv:2607.23074 [INSPIRE].

Notes

^1. Note that in contrast to ref. [1], we do not apply a long sequence of cooling maps with randomly chosen cycle times for each frequency, but instead vary both the cycle time and the bath frequency randomly.

^2. Recall that in our model, the energy of the ground state is always negative and the energy of the maximally mixed state is zero. Hence, the maximally mixed state in eq. (3.12) does not contribute to the energy and only a damping term appears.

^3. As mentioned before the analytical results derived in appendix A hold for any distribution. However, for the numerical derivations we choose here a discrete number of frequencies.

^4. Note that the result is the same as if we do it the other way around, see appendix A.1.

^5. Recall that the evolution of the mode pairs $(k,-k)$ decouples and, therefore, we can consider them separately.

^6. The reason for that is that both coefficients in front of $d_{\ell }$ and $q_{\ell }$ are smaller than 1 (for all but possibly a zero measure set of choices of the pairs $(\tau ,\mathrm{\Delta })$). Hence, in the limit $\ell \; \to \; \infty$ both, $d_{\ell }$ and $q_{\ell }$, will vanish. Setting $d_{\ell }=\; 0$ in the equation for $p_{\ell +1}$ we also have that the absolute value of $p$ decreases to $0$ and therefore $p_{\ell }$ vanishes too.

^7. This follows trivially from $\overline{n_{\ell }^2}-{\overline{n_{\ell }}}^2\le \overline{n_{\ell }^2}\le \overline{n_{\ell }}$.

^8. This state is then used as the starting state for adiabatic evolution, as detailed in appendix B.4.