# Quasiparticle cooling algorithms for quantum many-body state preparation

Jerome Lloyd  Department of Theoretical Physics, University of Geneva, Geneva,
Switzerland  Alexios Michailidis  Department of Theoretical Physics,
University of Geneva, Geneva, Switzerland  PlanQC GmbH, Lichtenbergstr. 8,
85748 Garching, Germany  Xiao Mi  Google Quantum AI, Santa Barbara CA, USA
Vadim Smelyanskiy  Google Quantum AI, Santa Barbara CA, USA  Dmitry A. Abanin
Google Quantum AI, Santa Barbara CA, USA  Department of Physics, Princeton
University, Princeton NJ 08544, USA

(28th August 2024)

###### Abstract

Probing correlated states of many-body systems is one of the central tasks for
quantum simulators and processors. A promising approach to state preparation
is to realize desired correlated states as steady states of engineered
dissipative evolution. A recent experiment with a Google superconducting
quantum processor [X. Mi et al., Science 383, 1332 (2024)] demonstrated a
cooling algorithm utilizing auxiliary degrees of freedom that are periodically
reset to remove quasiparticles from the system, thereby driving it towards the
ground state. We develop a kinetic theory framework to describe quasiparticle
cooling dynamics, and employ it to compare the efficiency of different cooling
algorithms. In particular, we introduce a protocol where coupling to
auxiliaries is modulated in time to minimize heating processes, and
demonstrate that it allows a high-fidelity preparation of ground states in
different quantum phases. We verify the validity of the kinetic theory
description by an extensive comparison with numerical simulations of a 1d
transverse-field Ising model using a solvable model and tensor-network
techniques. Further, the effect of noise, which limits efficiency of
variational quantum algorithms in near-term quantum processors, can be
naturally described within the kinetic theory. We investigate the steady state
quasiparticle population as a function of noise strength, and establish
maximum noise values for achieving high-fidelity ground states. This work
establishes quasiparticle cooling algorithms as a practical, robust method for
many-body state preparation on near-term quantum processors.

##  I Introduction

One of the central challenges for quantum simulation is to find reliable
algorithms to prepare desired entangled quantum many-body states. In
particular, preparing low-energy states of a given many-body Hamiltonian is an
essential task for probing quantum dynamics, modeling correlated quantum
materials, and for quantum chemistry applications of quantum computers [1, 2,
3, 4]. Classical tasks such as high-dimensional optimization problems can also
be formulated in terms of finding ground states of appropriate Hamiltonians
[5, 6]. A number of quantum algorithms for preparing ground states of many-
body Hamiltonians have been proposed [7, 8, 9, 10, 11, 12].

Two important classes of low-energy state-preparation protocols employed in
experiments with noisy intermediate-scale quantum (NISQ) devices [13] are: (i)
adiabatic protocols [9, 14], which rely on the adiabatic theorem and a slow
variation of Hamiltonian parameters, and have been widely used in analogue
quantum simulators such as ultracold atoms [15, 16] and Rydberg atom arrays
[17, 18]; (ii) variational methods [19, 20, 21], suitable for digital quantum
simulators, in which the target state is prepared by applying a variational
quantum circuit to a simple initial state (e.g. a product state).

Figure 1: (a) Schematic mechanism of quasiparticle cooling: quasiparticles in
the system, with energy above the system gap $\Delta$, are resonantly
transferred to auxiliary qubits, which are then reset to ground state,
removing entropy and energy from the system. (b) Digital quantum circuit used
in reset protocol: a system of $N_{S}$ qubits $\mathcal{Q}_{i}$ is coupled to
an auxiliary bath of $N_{A}$ qubits $\mathcal{A}_{i}$. Unitary evolution under
system, auxiliary, and interaction gates is applied, and repeated $T$ times,
followed by the reset of auxiliaries. This cooling cycle is repeated until the
system converges to its steady state. (c) Edge cooling setup. (d) Bulk cooling
setup, $N_{A}=N_{S}$.

The unitary state-preparation protocols on NISQ processors are hindered by the
presence of intrinsic decoherence. Such protocols have no native mechanism to
remove excitations, so decoherence will eventually thermalize the system to a
trivial state, e.g. a fully mixed state. This restricts the size of feasible
circuits and limits the duration for which a quantum state can be stored on
the device [22, 23]. Further, using variational circuits to prepare complex
quantum states may be challenging due to flat optimization landscapes (‘barren
plateaus’) [24, 25, 26].

A fundamentally different approach is to prepare many-body states as steady
states of engineered dissipative evolution [27]. Because the dissipation
directly competes with the decoherence effects, such state preparation methods
can partially overcome the difficulties posed by noise. However, understanding
which class of experimentally realisable dissipators can give rise to
interesting quantum correlated steady states is a highly non-trivial question
that has attracted much historic and recent interest [28, 29, 30, 31, 32, 33,
34]. In particular, specific protocols that involve local engineered
dissipation for frustration-free Hamiltonians (a restricted but important
class of models) have been proposed and implemented experimentally in few-
qubit systems of trapped ions and superconducting qubits [35, 36, 37].

More recently, in an experiment [38] on the digital Google Quantum Processor,
a physics-inspired cooling algorithm was demonstrated for cooling 1d and 2d
quantum spin systems. In this experiment, the system evolved under a periodic
(Floquet) quantum circuit specified by a unitary Floquet operator
$\hat{U}_{S}$, constructed as a shallow quantum circuit of gates with finite,
but small rotation angles. In this regime, the time evolution of local
observables is well-approximated by a prethermal Hamiltonian $H_{\rm eff}$
[39]. The approximate conservation of $H_{\rm eff}$ stems from the suppression
of heating processes, which are exponentially slow in the inverse rotation
angles of gates in $\hat{U}_{S}$. The key idea of the implemented algorithm
(displayed schematically in Fig. 1a,b — see also [40, 41, 42, 43, 44] for
related works) was to use a subset of qubits as auxiliaries, tuning their
quasienergies and couplings to system qubits in order to resonantly transfer
excitations from the system. The entropy of the system is removed by
periodically resetting the auxiliaries. In the experiment, it was shown that
the system reached a steady low-energy correlated state. Furthermore, the
achieved steady states could be efficiently purified to yield high fidelities
with respect to the true ground states, and it was argued that,
asymptotically, the dissipative preparation is advantageous compared to
variational unitary protocols.

This progress calls for the development of practical, robust, dissipative
many-body state preparation algorithms. In this work, we develop a theory of
quasiparticle cooling on quantum processors, concentrating on efficient
preparation of phases of matter with low-energy quasiparticle excitations, in
realistic noisy environments. In such many-body systems, excited states at
sufficiently low energy density can be described in terms of occupation
numbers of non-interacting bosonic or fermionic degrees of freedom [45]. The
quasiparticles generally have a finite lifetime, due to the decay to other
states; however, provided this lifetime is sufficiently long, quasiparticle
description can be employed. A celebrated example is the Landau Fermi-liquid
theory, which provides a description of interacting electrons in metals in
terms of fermionic quasiparticles [45]. We note that quasiparticle description
applies to a variety of phases of matter including Fermi liquids,
superconductors, and magnetically ordered states.

Focusing on 1d qubit – or, equivalently, spin-$1/2$ – systems, we consider
generalised auxiliary-reset algorithms encompassing the one implemented in the
experiment of Ref. [38]. We derive a dissipative kinetic equation describing
the quasiparticle population dynamics under the action of the protocol, which
allows for a transparent physical picture on how different parameters affect
the state preparation fidelity. In particular, we investigate the effects of:
(i) auxiliary layout as shown in Fig. 1c,d, with auxiliaries only at the
system boundaries (edge cooling setup) or with a finite density of auxiliaries
in the thermodynamic limit (bulk cooling setup); and (ii) the time dependence
of the system-auxiliary couplings. We find that the latter is crucial in
reaching high fidelities, and we introduce a protocol with modulated couplings
to drastically reduce unwanted heating processes arising due to the finite
duration of a cooling cycle. Below, we characterize the protocol performance
by log-fidelity per qubit, $-\log{\mathcal{F}}/N_{S}$, where $\mathcal{F}$ is
the fidelity of the steady state with respect to the true ground state, and
$N_{S}$ is the number of system qubits. We argue that, if the steady state is
sufficiently close to the ground sate, this quantity is related to the steady-
state quasiparticle density $n$ via $-\log{\mathcal{F}}/N_{S}\approx n$.
Within our kinetic equation approach, we find $n\sim\exp(-\Omega(\Delta T))$
under the modulated coupling protocol, where $\Delta$ is the quasiparticle
excitation gap and $T$ is the number of unitary layers before the auxiliaries
are reset (see Fig. 1b). Thus, high fidelity can be achieved by increasing
parameter $T$.

We assess the validity of our theory first with the example of an integrable
1d transverse-field Ising model (TFIM) spin chain, where quasiparticles are
infinitely long-lived and the relevant operators are known explicitly [46,
47]. This model exhibits two phases, a paramagnetic (PM) and an
antiferromagnetic (AFM) phase. The distinct nature of the quasiparticle type
in the two phases (magnons in the paramagnet and domain walls in the
antiferromagnet), as well as spontaneous symmetry breaking in the AFM phase,
lead to qualitative differences in the importance of different cooling
processes. By performing extensive numerical simulations, we show that the
kinetic equation correctly captures the system’s evolution and its steady
state even at a relatively strong auxiliary-system coupling — the regime
relevant for current experiments. We compare the experimental protocol of Ref.
[38] (where the couplings were kept constant) to our modulated coupling
protocol, finding that the latter allows for much higher fidelities in the
limit of large $T$.

Noise, inevitably present in today’s quantum simulators and processors,
severely limits the performance of quantum algorithms. For example,
superconducting qubits are subject to dephasing and decay processes, which
introduce unwanted quasiparticles to the system. We modify the kinetic
equation by incorporating noise-induced heating on top of the cooling
protocol, and investigate how the non-equilibrium distribution of
quasiparticle level occupations in the steady state is affected by noise. For
the Floquet TFIM example, we observe that the density $n$ of quasiparticles in
the steady state scales differently with the Markovian noise strength in the
two phases. Characterizing decoherennce by a dimensionless parameter $\gamma$,
which reflects the noise-induced error probability per qubit over one period
of Floquet evolution, we find scaling $n\propto\gamma$ in the PM phase, while
in the AFM phase, we observe a crossover between linear behaviour and a
scaling $n\propto\sqrt{\gamma}$. This scaling behaviour was predicted in Ref.
[43], and can be understood within the kinetic framework from the nature of
different quasiparticle types in the two phases. Achieving low quasiparticle
densities in the AFM phase therefore requires significantly weaker noise,
compared to the PM phase. Nevertheless, we find that the modulated coupling
protocol achieves high fidelities when the noise rate is small compared to the
auxiliary-system coupling strength.

While we largely focus on cooling in the integrable TFIM model, we expect the
quasiparticle cooling algorithm to be effective also for non-integrable
systems. Indeed, many such systems have low-energy quasiparticles which become
weakly interacting at low energy densities. Thus, at least at low energies
(the most relevant regime for ground state preparation), our analysis should
apply. We illustrate the versatility of our approach by studying cooling in a
TFIM with broken integrability [48], where we observe similar efficiency of
our protocol as in the integrable case.

The rest of the paper is organized as follows: in Section II we describe the
quasiparticle cooling protocol, and derive an equation that describes the
evolution of the quasiparticle level populations, assuming a weak auxiliary-
system coupling. We introduce the modulated coupling protocol (MCP) to
suppress unwanted heating processes. In Section III, we turn to cooling the
integrable TFIM. We analyse cooling with boundary auxiliaries and with a
finite density of auxiliaries, for different cooling protocols. We find that
the MCP leads to effective cooling, which can be accurately captured by the
kinetic equation of Sec. II. Noise is introduced in Section IV: we examine how
noise processes (dephasing and decay) affect the cooling, and show that these
effects can be included into the kinetic description without complication. In
Section V, we apply our protocol to study a non-integrable model, the TFIM
with an additional longitudinal field. Finally, in Section VI, we provide
concluding remarks and directions for future work.

##  II Quasiparticle cooling protocol

###  II.1 Setup and protocol

We consider the following general setup for a system of qubits, or
spins-$1/2$, a particular version of which was realized in the recent
experiment with the Google quantum processor [38]. The central idea follows
Ref. [27], and the main developments we present here are (i) the derivation of
a simplified kinetic theory, based on the quasiparticle picture, and (ii) the
introduction of time-dependent system-bath coupling, which, as we will see
below, is crucial to reach high fidelities.

The qubits are divided into two groups, $N_{S}$ qubits forming a system of
interest, and $N_{A}$ auxiliary qubits acting as a bath. One cycle of the
protocol (_cooling cycle_), which effectively realizes a Floquet system
coupled to an engineered bath, entails the following steps:

  * (i)

Auxiliary qubits are initialized in a state described by a density matrix
$\hat{\rho}_{B}^{0}$.

  * (ii)

$T$ layers of unitary evolution are applied. A unitary applied at time step
$1\leq\tau\leq T$ is a combination of the Floquet operator $\hat{U}_{S}$
acting on system qubits, a unitary $\hat{U}_{B}$ acting on the bath qubits,
and an interaction term $\hat{U}_{\theta,\tau}$ which couples system and bath
qubits. The interaction terms is generally time-dependent, and $\theta$ is the
coupling strength (see Eq. (8) below).

  * (iii)

The bath qubits are reset into the initial state $\hat{\rho}_{B}^{0}$, which
we will assume to be a stationary state of the bath unitary $\hat{U}_{B}$.

The main protocol control parameters are $T$, which we refer to as the _reset
time_ , and the system-auxiliary coupling. Resetting auxiliary qubits at the
end of the cycle defines a quantum channel acting on the system,

$\Phi(\hat{\rho}_{S})=\text{Tr}_{B}\
\hat{U}_{T}(\hat{\rho}_{S}\otimes\hat{\rho}_{B}^{0})\hat{U}_{T}^{\dagger},$
(1)

where $\hat{U}_{T}$ is the combined cooling cycle unitary,

$\hat{U}_{T}=\hat{U}_{\theta,T}\hat{U}_{B}\hat{U}_{S}\ldots\hat{U}_{\theta,1}\hat{U}_{B}\hat{U}_{S},$
(2)

the operator $\hat{\rho}_{S}$ is the system density matrix at the start of the
cycle, and $\text{Tr}_{B}$ is the trace over the bath degrees of freedom. The
system density matrix after $\nu$ cooling cycles is given by

$\hat{\rho}_{S}^{\nu}=\Phi^{\nu}(\hat{\rho}_{S}^{0}).$ (3)

We will assume the steady state of the dynamics is unique, and is therefore
reached independent of the initial state. In the analysis below we therefore
take the initial state to be the maximally mixed state,
$\hat{\rho}_{S}^{0}\propto\hat{I}$, unless specified otherwise.

The system Floquet $\hat{U}_{S}$ can in particular correspond to an exact or
prethermal effective Hamiltonian $\hat{H}_{\rm eff}$,

$\hat{U}_{S}\approx\exp\bigg{(}-\frac{i\pi}{2}\hat{H}_{\rm eff}\bigg{)},$ (4)

and our goal is to engineer dissipative evolution such that the steady state
($\nu\to\infty$) of the system is close to the ground state of $\hat{H}_{\rm
eff}$. In order to have a concrete model in mind, we will later study cooling
in the case of the 1d Floquet transverse-field Ising model,

$\displaystyle\hat{U}_{S}=\exp\bigg{(}{-\frac{i\pi
J}{2}\sum_{i=1}^{N_{S}-1}\hat{X}_{i}\hat{X}_{i+1}}\bigg{)}\exp\bigg{(}\frac{i\pi
g}{2}\sum_{i=1}^{N_{S}}\hat{Z}_{i}\bigg{)}.$ (5)

Here and below $\hat{X}_{i},\hat{Y}_{i},\hat{Z}_{i}$, denote the Pauli
operators for qubit $i$, normalized as
$\hat{X}_{i}^{2}=\hat{Y}_{i}^{2}=\hat{Z}_{i}^{2}=1$. The Floquet TFIM is
exactly solvable via the Jordan-Wigner mapping to a model of non-interacting
fermionic quasiparticles (see Appendix C), and so provides a convenient
setting to study quasiparticle cooling.

We keep the bath evolution simple, assuming for the initial bath state the
product state with all spins up,

$\hat{\rho}_{B}^{0}=\otimes_{j=1}^{N_{A}}|0\rangle_{j}\,{}_{j}\langle 0|$ (6)

— note that in our convention the spin-up state is identified with the
$\ket{0}$ state of the qubits — and for the bath unitary we choose

$\hat{U}_{B}=\exp\bigg{(}\frac{i\pi
h}{2}\sum_{j=1}^{N_{A}}\hat{Z}_{j}\bigg{)},$ (7)

so that each auxiliary qubit simply acquires a phase in the eigenbasis of the
$\hat{Z}$ operator. More complicated (e.g. interacting) baths can be
considered with possible advantages. Parameter $\pi h$ is referred to as the
_auxiliary quasienergy_ and does not depend on the unitary step $\tau$.

In contrast, the system-bath coupling will in general vary in time:

$\hat{U}_{\theta,\tau}=\exp(-i\theta f_{\tau}\hat{V}),\hskip
14.22636pt\hat{V}=\sum_{j=1}^{N_{A}}\hat{V}_{j}.$ (8)

We have chosen to separate the time-dependence into the function $f_{\tau}$,
normalised according to $\sum_{\tau=1}^{T}f_{\tau}=1$ 111For protocols
considered below, $f_{\tau}$ can be positive or negative, but for this choice
the sum $\sum_{\tau=1}^{T}|f_{\tau}|\gtrsim 1$., keeping $\theta$ as a small
parameter which controls the coupling strength. Both $f_{\tau}$ and $\theta$
are dimensionless, and the interaction operator, which has dimension of
energy, can in general be written in the form

$\hat{V}=\sum_{a}\hat{A}^{a\dagger}\hat{B}^{a}+\hat{A}^{a}\hat{B}^{a\dagger},$
(9)

where operators $\hat{A}$ ($\hat{B}$) acts on system (bath). Note that index
$a$ does not have to be identical to the index numbering auxiliary qubits.

A particular form of interaction that we will consider below arises when an
auxiliary $\mathcal{A}_{j}$ is coupled to a single qubit $\mathcal{Q}_{j}$ of
the system by an operator $\hat{V}$ that can generate a partial iSWAP gate:

$\hat{V}=\pi\sum_{j}(\hat{\sigma}^{+}_{\mathcal{Q}_{j}}\hat{\sigma}^{-}_{\mathcal{A}_{j}}+\hat{\sigma}^{-}_{\mathcal{Q}_{j}}\hat{\sigma}^{+}_{\mathcal{A}_{j}}),$
(10)

where $\mathcal{Q}_{j}$ ($\mathcal{A}_{j}$) represent system (auxiliary)
position associated with $j$-th coupling. Here
$\hat{\sigma}^{\pm}=\frac{1}{2}(\hat{X}\pm i\hat{Y})$ are spin
raising/lowering operators.

###  II.2 Weak-coupling analysis

In the remaining theoretical analysis, we assume that the excitations of the
many-body system are well-defined, long-lived quasiparticles, and that we are
in the limit of weak system-bath coupling, $\theta|V_{j}|\ll 1$. The
quasiparticle assumption holds exactly in the example of the TFIM, Eq. (5); we
note, however, that interacting systems away from criticality typically have
quasiparticle excitations – examples include Fermi-liquids, superconductors,
and symmetry-broken magnetic phases [45]. We expect the quasiparticle cooling
algorithm to remain effective for such interacting systems. Our strategy will
be to derive an equation describing evolution of the quasiparticle level
occupations in the weak-coupling limit.

We define the quasiparticle occupation number operator $\hat{n}_{k}$, where
index $k$ labels fermionic quasiparticle levels (e.g. momentum in a system
with periodic boundary conditions). The quasiparticles are sufficiently long-
lived, such that $\hat{n}_{k}$ is an integral of motion of the system’s
evolution operator,

$\hat{n}_{k}\hat{U}_{S}=\hat{U}_{S}\hat{n}_{k}.$ (11)

In practice, it is sufficient to assume that the quasiparticle lifetime is
much longer than the duration of one cooling cycle.

It is convenient to work in the interaction picture with respect to the
evolution operator $\hat{U}_{0}=\hat{U}_{B}\hat{U}_{S}$, such that an operator
$\hat{O}_{\tau}$ (originally in the Schrödinger picture) is given by:

$\hat{O}_{I\tau}=\hat{U}_{0}^{-\tau}\hat{O}_{\tau}\hat{U}_{0}^{\tau},\;\;\;1\leq\tau\leq
T.$ (12)

The evolution operator can be rewritten as follows,

$\hat{U}_{T}=\hat{U}_{0}^{T}\mathcal{\hat{U}}_{T},\;\;\;\mathcal{\hat{U}}_{T}=\mathcal{T}\exp\Big{(}-i\theta\sum_{\tau=1}^{T}f_{\tau}\hat{V}_{I\tau}\Big{)},$
(13)

where $\mathcal{T}$ denotes time-ordering. Taking into account Eqs. (1,11),
the quasiparticle occupation number after one cooling cycle reads:

$n_{k}^{\prime}={\rm
Tr}_{S}(\Phi(\hat{\rho}_{S})\hat{n}_{k})=\text{Tr}_{SB}\Big{(}\hat{\rho}_{S}\otimes\hat{\rho}_{B}^{0}\Big{[}\mathcal{\hat{U}}^{\dagger}_{T}\hat{n}_{k}\mathcal{\hat{U}}_{T}\Big{]}\Big{)},$
(14)

At weak coupling, the occupation number change over one cycle, $\delta
n_{k}=n_{k}^{\prime}-n_{k}$, is small, and we compute it by expanding the
r.-h.s. of the above equation to the second order in $\theta$. Noting that the
linear order in $\theta$ vanishes, and defining

$\delta n_{k}={\rm Tr}_{S}(\hat{\rho}_{S}\delta\hat{n}_{k}),\hskip
14.22636pt\delta\hat{n}_{k}=\text{Tr}_{B}\Big{(}\hat{\rho}_{B}^{0}\Big{[}\mathcal{\hat{U}}^{\dagger}_{T}\hat{n}_{k}\mathcal{\hat{U}}_{T}\Big{]}\Big{)}-\hat{n}_{k},$
(15)

this yields:

$\displaystyle\delta\hat{n}_{k}=\theta^{2}\sum_{ab}\sum_{\tau_{1}=1}^{T}\sum_{\tau_{2}=1}^{T}f_{\tau_{1}}f_{\tau_{2}}\Gamma_{ab}(\tau_{1}-\tau_{2})\Big{(}\hat{A}^{a\dagger}_{I\tau_{1}}\hat{n}_{k}\hat{A}^{b}_{I\tau_{2}}-\frac{1}{2}\\{\hat{n}_{k},\
\hat{A}^{a\dagger}_{I\tau_{1}}\hat{A}^{b}_{I\tau_{2}}\\}\Big{)}+i\theta^{2}[\hat{H}_{\rm
LS},\hat{n}_{k}].$ (16)

Here we defined the autocorrelators of the $\hat{B}$ operators,

$\Gamma_{ab}(\tau)=\text{Tr}_{B}\Big{(}\hat{\rho}_{B}^{0}\hat{B}^{a\dagger}_{I\tau}\hat{B}^{b}_{I0}\Big{)},$
(17)

and used relations $\Gamma_{ab}(\tau)=\Gamma_{ba}^{*}(-\tau)$. The second term
in the r.-h.s. of Eq. (16) can be viewed as arising from the “Lamb shift”
renormalization:

$\displaystyle\hat{H}_{\rm
LS}=\frac{i}{2}\sum_{ab}\sum_{\tau_{1}=1}^{T}\sum_{\tau_{2}=1}^{T}f_{\tau_{1}}f_{\tau_{2}}(\Theta(\tau_{1}-\tau_{2})-\Theta(\tau_{2}-\tau_{1}))\Gamma_{ab}(\tau_{1}-\tau_{2})\hat{A}^{a\dagger}_{I\tau_{1}}\hat{A}^{b}_{I\tau_{2}}.$
(18)

where $\Theta(\tau)$ is the standard Heaviside step function, with the
convention $\Theta(0)=1/2$.

###  II.3 Rate equation

Note that the above derivation so far follows the steps of the textbook
derivation of the Lindblad equation [50, 51], except here we focus on the
Heisenberg-picture evolution of the occupation number operators, rather than
on the system’s density matrix. However, to derive the Lindblad equation one
typically assumes a rapid decay of bath temporal correlations, compared to
system’s dynamical timescales; in contrast, in our analysis the memory effects
of the bath will play an essential role.

To simplify the analysis, we make an assumption that the system’s density
matrix in the beginning of a cooling cycle is diagonal in the quasiparticle
basis, $[\hat{\rho}_{S},\hat{n}_{k}]=0$. Moreover, we will make a further
assumption that the correlations between occupations of different levels can
be neglected (closely related to molecular chaos assumption in kinetic
theory). The former assumption can be justified in the limit $\theta\to 0$, by
coarse-graining over sufficiently many cooling cycles (secular approximation).
The latter assumption is more difficult to prove 222We note that the second
assumption is practically necessary, as the diagonal assumption still leaves
us with $O(\exp(N_{S}))$ coupled equations to solve.. However, as we will show
below, these assumptions lead to a theory of quasiparticle cooling which is in
quantiative agreement with numerical simulations.

Under these two assumptions, the system’s density matrix is given by
$\hat{\rho}_{S}=\sum_{\vec{\alpha}}\rho_{\vec{\alpha}\vec{\alpha}}\ket{\vec{\alpha}}\bra{\vec{\alpha}}$,
where $\ket{\vec{\alpha}}$ form a basis of states with occupation numbers
$\\{\alpha_{k}\\}$, $\alpha_{k}=0,1$ and quasienergy

$\varepsilon(\vec{\alpha})=\sum_{k}\epsilon_{k}\alpha_{k},$

where $\epsilon_{k}$ is quasienergy of quasiparticle level $k$. The
probability $\rho_{\vec{\alpha}\vec{\alpha}}$ can be expressed in terms of
average occupation numbers $n_{q}={\rm Tr}(\hat{\rho}_{S}\hat{n}_{q})$ via

$\rho_{\vec{\alpha}\vec{\alpha}}=\prod_{q}n_{q}^{\alpha_{q}}(1-n_{q})^{1-\alpha_{q}}.$
(19)

Note that in writing this form, we used the assumption that there are no
correlations between occupations of different quasiparticle states. This
ansatz can also be viewed as a time-dependent Generalized Gibbs ensemble (GGE)
ansatz for $\hat{\rho}$ [53, 54].

Next, we compute the change in the occupation number over one cycle, $\delta
n_{k}$ using the above form of the density matrix and Eq. (16). Since density
$\hat{n}_{k}$ commutes with $\hat{\rho}_{S}$, the Lamb shift contribution
plays no role. Denoting the Bohr frequencies by

$\Delta(\vec{\alpha},\vec{\beta})=\varepsilon({\vec{\alpha}})-\varepsilon({\vec{\beta}}),$
(20)

we obtain

$\delta
n_{k}=\theta^{2}\sum_{\vec{\alpha},\vec{\beta}}\sum_{a,b,\tau_{1},\tau_{2}}\rho_{\vec{\alpha}\vec{\alpha}}(\beta_{k}-\alpha_{k})f_{\tau_{1}}f_{\tau_{2}}\Gamma_{ab}(\tau_{1}-\tau_{2})e^{i(\tau_{1}-\tau_{2})\Delta(\vec{\alpha},\vec{\beta})}\bra{\vec{\alpha}}\hat{A}^{a\dagger}\ket{\vec{\beta}}\bra{\vec{\beta}}\hat{A}^{b}\ket{\vec{\alpha}}$
(21)

This equation illustrates that the quasiparticle number evolution is
controlled by (i) spectral correlations of the bath, given by
$\Gamma_{ab}(\tau_{1}-\tau_{2})$, (ii) the choice of the time-dependent
coupling, $f_{\tau}$, and (iii) the matrix elements of the operators
$\hat{A}^{a}$ entering auxiliary-system coupling. With our above choice of
bath unitary (7) and coupling (10), the bath spectral function reads

$\Gamma_{\mathcal{A}_{j}\mathcal{A}_{l}}(\tau)=\bra{\vec{0}}\hat{\sigma}^{-}_{\mathcal{A}_{j}}(\tau)\hat{\sigma}^{+}_{\mathcal{A}_{l}}(0)\ket{\vec{0}}=\delta_{jl}e^{-i\pi
h\tau},$ (22)

and Eq. (21) takes the following form:

$\delta
n_{k}=\theta^{2}\sum_{j=1}^{N_{A}}\sum_{\vec{\alpha}\vec{\beta}}\rho_{\vec{\alpha}\vec{\alpha}}(\beta_{k}-\alpha_{k})|F_{h,T}(\Delta(\vec{\alpha},\vec{\beta}))|^{2}|\bra{\vec{\beta}}\hat{\sigma}_{\mathcal{Q}_{j}}^{+}\ket{\vec{\alpha}}|^{2}.$
(23)

where

$F_{h,T}(\epsilon)=\pi\sum_{\tau=1}^{T}f_{\tau}e^{i\tau(\epsilon-\pi h)}.$
(24)

We refer to $F_{h,T}(\epsilon)$ as the _filter function_ associated with
$f_{\tau}$, essentially being its Fourier transform. Eq. (23) is a general
rate equation for quasiparticle level occupation numbers. Next, we discuss how
the auxiliary quasienergy $\pi h$ and the coupling $f_{\tau}$ can be tailored
to remove quasiparticles from the system, driving it towards the ground state.

###  II.4 Two cooling protocols

In Eq. (23) above we expressed the change in the occupation numbers in terms
of the transition rates between many-body states $|\vec{\alpha}\rangle$ and
$|\vec{\beta}\rangle$ induced by the coupling to auxiliaries. It is natural to
consider the prethermal limit where quasiparticle quasienergies are close to
the energies of quasiparticles of the corresponding prethermal Hamiltonian
[39]. For example, transverse-field Ising model with arbitrary integrability-
breaking local terms is in the prethermal regime as long as $|J|,|g|\ll 1$ and
a similar condition is satisfied by the additional, integrability-breaking
terms. In the prethermal regime, generally $|\epsilon_{k}|\ll\pi$. We note
that in this limit, the heating processes due to Floquet evolution, and the
associated energy non-conservation, are suppressed exponentially in
$1/|g|,1/|J|$.

The above rate equation (23) includes both heating
($\Delta({\vec{\alpha},\vec{\beta}})<0$) and cooling
($\Delta({\vec{\alpha},\vec{\beta}})>0$) processes. To ensure effective
cooling, the filter function should be engineered to satisfy (for
$\epsilon>0$):

$\bigg{|}\frac{F_{h,T}(+\epsilon)}{F_{h,T}(-\epsilon)}\bigg{|}\gg 1.$ (25)

The prethermalization assumption allows us to neglect processes where
quasienergy changes by $\pm 2\pi,\pm 4\pi$ etc, since the corresponding matrix
elements are strongly suppressed. Below, we will introduce two protocols that
realise Eq. (25) to different degrees.

####  II.4.1 Step-wise cooling protocol (SCP)

In the first strategy, one aims to resonantly remove single quasiparticles
from the system, by approximately matching auxiliary excitation quasienergy
$\pi h$ to a value $\epsilon$ within the quasiparticle band,
$\Delta<\epsilon<\Delta+W$, where $\Delta$ is the single-particle quasienergy
gap and $W$ is the bandwidth. The function $f_{\tau}=1/T$ (26) is kept
constant, therefore giving a ‘step’ pulse for $1\leq\tau\leq T$. We will refer
to this approach as the step-wise cooling protocol (SCP).

We note that the SCP is the strategy used in Ref. [38]. The filter function
corresponding to the SCP gives, up to a phase factor,

$F^{\rm SCP}_{h,T}(\epsilon)=\frac{\pi}{T}\frac{\sin[(\epsilon-\pi
h)T/2]}{\sin[(\epsilon-\pi h)/2]},$ (27)

which for $T\gg 1$ has a sharp peak at $\epsilon=\pi h$, with a width of order
$\delta\epsilon=\frac{2\pi}{T}$. Thus, in the SCP, the processes of
quasiparticle removal with quasienergies $\epsilon\approx\pi h$ are favored
approximately according to

$\bigg{|}\frac{F^{\rm SCP}_{h,T}(+\pi h)}{F^{\rm SCP}_{h,T}(-\pi
h)}\bigg{|}\approx hT.$ (28)

The SCP is efficient in systems with a sufficiently large quasiparticle gap,
as demonstrated in [38]. However, as the gap is decreased, in particular when
a quantum critical point is approached, the heating processes become more
important. Even in a gapped phase, the suppression of heating is only
algebraic in quasiparticle energy, which leads to non-zero quasiparticle
density in the steady state. Therefore, SCP does not allow to prepare a many-
body ground state with high fidelity.

####  II.4.2 Modulated cooling protocol (MCP)

Alternatively, the coupling strength $f_{\tau}$ can be modulated in time such
that $F_{h,T}(\epsilon)$ leads to approximately exponential suppression of
heating transitions. One particular choice of $f_{\tau}$ is
$f_{\tau}=\frac{1}{A\beta}\frac{\sin[\pi(\tau-\tau_{0})/2]}{\sinh[\pi(\tau-\tau_{0})/\beta]},$
(29) where $\tau_{0}=T/2$, $A$ is the normalization factor 333We note that the
normalization factor is chosen according to the condition
$\sum_{\tau=1}^{T}f_{\tau}=1$. The modulated coupling in Eq. (29) is sign-
changing, however, $\sum_{\tau=1}^{T}|f_{\tau}|$ is of order 1 for all cases
considered below.. The auxiliary quasienergy is fixed at $\pi h=\pi/2$. We
call this approach the modulated cooling protocol (MCP).

In the SCP, transitions occur at the natural resonant frequency of the
auxiliary qubits; in the MCP, the resonant frequencies are engineered
according to modulation in the time-domain. More specifically, the MCP is
designed to satisfy the thermal detailed balance condition, and $F_{h,T}$
takes the form of a broadened step-function (Fermi distribution) in frequency
space for $1\ll\beta\ll T$, with ‘thermal’ broadening $\sim 1/\beta$:

$F^{\rm
MCP}_{h,T}\approx\pi\bigg{[}\frac{\tanh(\epsilon\beta/2)-\tanh((\epsilon-\pi)\beta/2)}{\tanh(\pi\beta/4)}\bigg{]}.$
(30)

We refer the reader to Appendix A for the derivation of this formula in the
limit where $T\to\infty$, and $0<\beta/T\ll 1$, and Appendix B for a more
detailed discussion of the detailed balance condition. The frequency-step form
of the MCP leads to a suppression

$\bigg{|}\frac{F^{\rm MCP}_{h,T}(+\epsilon)}{F^{\rm
MCP}_{h,T}(-\epsilon)}\bigg{|}\approx e^{\epsilon\beta},$ (31)

and heating transitions are now _exponentially_ suppressed. We have also
verified that the step-like form of $F^{\rm MCP}_{h,T}(\epsilon)$ is realized
for moderate values of $T\approx 10$, which is the practical regime in the
presence of noise.

In the rest of this paper, we focus on cooling _gapped_ quasiparticle systems,
with $\Delta>0$. As in the case of thermal equilibrium, we anticipate the
difficuly of cooling gapless systems to their ground state, since zero energy
modes may be excited by arbitrarily small heating elements for any finite
$\beta$ (in the case of the MCP). Indeed, on general grounds of the energy-
time uncertainty relation, we expect modes of energy $\epsilon$ require times
$T\approx\epsilon^{-1}$ to cool.

We note that other, more general protocols can be constructed, which in
certain cases may improve cooling efficiency. In particular, $h$ can be
modulated during the cooling cycle, as in the adiabatic demagnetization
protocol recently studied in Ref. [43]. In Sec. V, when considering cooling of
non-integrable systems, we will also vary the coupling strength $\theta$ over
cooling cycles in order to reach the steady state more rapidly. Our intention
with introducing the two protocols above is to provide a basis for effective
cooling, and in particular to highlight the role of modulating the coupling
pulse to block unwanted heating.

##  III Theory of quasiparticle cooling in a spin chain

In this Section, we study the performance of the cooling algorithms using the
example of the TFIM, specified by the Floquet operator in Eq. (5). We assume
the absence of noise (which we include in Section IV). We first consider the
edge auxiliary setup (Fig. 1c). We will see that this leads to solvable
dynamics in the sense that the density matrix of the system remains Gaussian.
We show that the predictions of the rate equation derived above are in
excellent agreement with the exact numerical solution. Further, we study the
bulk auxiliary setup (Fig. 1d), where integrability is broken by the coupling
to auxiliaries. In this case we compare the rate equation predictions to
matrix-product states simulations of the dissipative evolution, finding again
a good agreement. We investigate the dominant cooling and heating processes in
the two phases, paramagnetic (PM) and antiferromagnetic (AFM), of the TFIM,
finding qualitative differences due to the different nature of quasiparticles.
By comparing the performance of the SCP and MCP, we find that the latter
allows ground state preparation with a fidelity arbitrarily close to 1, in the
noiseless case, by increasing the parameter $\beta$.

###  III.1 Transverse-field Ising model

We start by briefly summarizing the properties of the Floquet TFIM. This model
can be mapped onto a $p$-wave superconducting chain by the Jordan-Wigner
transformation:

$\hat{\sigma}_{j}^{-}=\prod_{i=1}^{j-1}e^{i\pi\hat{c}_{i}^{\dagger}\hat{c}_{i}}\hat{c}_{j}^{\dagger},\hskip
14.22636pt\hat{Z}_{j}=1-2\hat{c}_{j}^{\dagger}\hat{c}_{j},$ (32)

where $\hat{c}_{j}^{\dagger},\hat{c}_{j}$ are fermionic creation/annihilation
operators. The quasiparticle bands have quasienergies $\epsilon_{k}$ given by:

$\cos\epsilon_{k}=\cos(\pi J)\cos(\pi g)-\sin(\pi J)\sin(\pi g)\cos k,$ (33)

where $k\in[-\pi;\pi]$ is a quasimomentum. As a function of parameters $J$ and
$g$, the model has four phases, separated by phase transition lines at $J=g$,
$J=1-g$, where the quasiparticle gap $\Delta$ closes. The four phases are
distinguished by the $0$ and $\pi$ Majorana edge modes [56, 57, 58]. Here we
will be interested in the part of the phase diagram where $0<J,g<1/2$, and two
phases exist: a paramagnetic (PM) phase at $g>J$, and an antiferromagnetic
(AFM) phase at $g<J$. For $g,J\to 0$ these are the phases of the Hamiltonian
TFIM; by analogy to this case, we define the Floquet vacuum projector

$\Omega=\sum_{g}|\Omega_{g}\rangle\langle{\Omega_{g}}|$ (34)

to project onto states where the valence quasiparticle band is filled, and the
conduction band is empty i.e. $\hat{c}_{k}|\Omega_{g}\rangle=0$, $\forall k$.
The index $g$ runs over the possible ground states: in the PM phase the vacuum
state is unique, while in the AFM phase in the limit $N_{S}\to\infty$, the
vacuum state is doubly degenerate, $g=\pm 1$, and the $\mathbb{Z}_{2}$
symmetry is spontaneously broken. In the fermionic counterpart of the TFIM,
this degeneracy stems from the existence of zero-(quasi)energy Majorana edge
modes, and the two vacuum states of the AFM differ by occupation of the
Majorana level. The quasienergy gap between the ground and first excited state
is given by $\Delta=2\pi|J-g|$. In the rest of the paper, unless otherwise
stated, results for the AFM (PM) phase correspond to the parameter points
$J=0.2\ (0.1),\ g=0.1\ (0.2)$.

For the open-boundary chain we consider, the allowed quasimomenta are
quantized, $k_{m}\approx\pi(m-1)/N_{S}$, $m=1,2,..,N_{S}$, and the fermionic
eigenmodes, defined by creation/annihilation operators
$\hat{\eta}_{k}^{\dagger},\hat{\eta}_{k}$, are superpositions of plane waves
with wave-vectors $\pm k$ (standing waves). The relation between on-site
fermionic operators and eigenmodes is

$c_{j}=\sum_{k}u_{jk}\eta_{k}+v_{jk}\eta_{k}^{\dagger},$ (35)

with the coefficients $u$, $v$ referred to as the ‘Bogolioubov coefficients’.
We refer to Appendix C for a more detailed discussion and explicit relations.

###  III.2 A solvable model: cooling at the edge

Figure 2: Performance of edge cooling protocol, for the parameter choice
described in the main text. (a) Evolution of quasiparticle density $n(t)$ as a
function of rescaled timescale $t=\nu\theta^{2}$, where $\nu$ is total number
of cooling cycles, for the PM phase. We compare exact numerical results with
$\theta=0.02$ (solid lines) and the rate equation (dashed lines), for both SCP
and MCP. (b) Same as in (a) but for AFM phase; (c) Log-fidelity per qubit of
the steady state with quasiparticle vacuum, obtain from kinetic equation, as a
function of the MCP parameter $T=3\beta$. Black squares are exact numerical
data, with $\theta=0.001$, showing that kinetic theory accurately captures the
steady-state properties.

We now study a solvable example of cooling in the edge auxiliary setup of Fig
1c. In this case, the auxiliary coupling operator $\hat{\sigma}_{1}^{\pm}$ is
free of the Jordan-Wigner string in the fermionic language, being proportional
to $c_{1}^{\dagger},c_{1}$ (see Eq. (32)). Therefore, in the rate equation
(23), only processes where occupation of one of the quasiparticle levels
changes by $\pm 1$ are allowed. Due to reflection symmetry, the auxiliary at
the other edge contributes equally to the cooling process, and after some
algebra, Eq. (23) takes the form:

$\displaystyle\delta n_{k}=-n_{k}W_{k}^{-}+(1-n_{k})W^{+}_{k},$ (36)

with the single-particle cooling and heating rates

$\displaystyle W^{-}_{k}=2\theta^{2}|F_{h,T}(+\epsilon_{k})|^{2}|u_{1k}|^{2},$
(37) $\displaystyle
W^{+}_{k}=2\theta^{2}|F_{h,T}(-\epsilon_{k})|^{2}|v_{1k}|^{2}.$ (38)

The coefficients $|u_{1k}|^{2},|v_{1k}|^{2}$ naturally represent the
probabilities for the quasiparticle (quasihole) to be found at the boundary.
Assuming that the system is away from the critical point, the probabilities
scale as $O(N_{S}^{-1})$ (from the normalization condition) away from the band
edge, but are $O(k^{2}N_{S}^{-1})\approx O(N_{S}^{-3})$ near the band edge
$k\approx 0,\pi$. This leads to a pronounced dependence of the cooling rates
on the quasienergy $\epsilon_{k}$, with the late time dynamics dominated by
the cooling of modes near the band edges. Eq. (36) can be solved for the
steady state populations

$n_{k}^{\infty}=\frac{W^{+}_{k}}{W^{+}_{k}+W^{-}_{k}},$ (39)

showing that a finite heating rate will inevitably lead to non-zero
quasiparticle population.

We will compare the rate equation predictions to exact numerical simulations,
which can be performed efficiently in the case of edge cooling: due to the
free-fermion nature of the TFIM with boundary auxiliaries, Eq. (1) defines a
Gaussian channel, so starting from an initial Gaussian state, such as the
identity state, Wick’s theorem guarantees that the complete many-body state
information is captured by the $2N\times 2N$-dimensional matrix of 2-point
fermion correlators, where $N=N_{S}+2$ is the number of system qubits plus the
two auxiliaries. By numerical evolution of the correlation matrix we can
assess the accuracy of the rate equation in capturing the two protocols.
Further detail on the numerical steps is provided in Appendix E.

As a first characterisation of the cooling performance, we study the dynamics
of the quasiparticle density $n(t)=N_{S}^{-1}\sum_{k}n_{k}(t)$. In Fig. 2a,b
we plot $n(t)$ for both protocols, SCP and MCP, in the two phases, starting
from the maximally mixed state $n(0)=1/2$. We take $T=4$, $h=0.3$ in the SCP,
and $T=28=3\beta$, $h=0.5$ in the MCP; the smaller reset time in the SCP case
is approximately optimal (see Fig. 2c) while the MCP choice corresponds to
$\beta\approx 2\pi\Delta^{-1}$. In both cases $\theta=0.02$ and $N_{S}=20$. We
observe a very good agreement between the rate equation theory (dashed lines)
and the exact numerical simulation (solid lines). In the considered time
window, for both PM and AFM phases under the MCP, the density decays
predominantly exponentially in time, with a rate proportional to $\theta^{2}$.
In accordance with the arguments regarding heating rates in Section II, the
MCP exhibits a saturation of $n(t)$ at a much smaller value than the SCP
(though in that case, the saturation is still to a relatively small value on
account of the large quasiparticle gap). We also note the long time-scales
associated with cooling to states of low quasiparticle density, with
$\nu\approx 300\theta^{-2}$ for the PM phase under the MCP and even greater
for the AFM. This is a consequence of the edge cooling setup, and the scaling
of the Bogoliubov coefficients near the band edge discussed above.

As a second test, we examine the fidelity of the steady state distribution
with respect to the Floquet vacuum, $\mathcal{F}=\text{Tr}(\Omega\rho).$ In
the case where Eq. (19) holds, the fidelity can be rewritten

$\mathcal{F}=\prod_{k}(1-n_{k})\approx\exp(-\sum_{k}n_{k}).$ (40)

The second line holds in the limit of sufficiently low total quasiparticle
density. In this regime, the log-fidelity per qubit and quasiparticle density
are approximately synonymous. In Fig. 2c we plot the steady state fidelity for
the two pulses as the parameter $\beta$ is increased, with $T=3\beta$ fixed
(for the SCP this is just a variation of parameter $T$). For parameters
considered, the SCP allows to reach fidelity of order $0.99$ per qubit while
the MCP fidelity (per qubit) is better than one part in $10^{6}$ by the point
$\beta\Delta\approx 4\pi$. From Eq. (31), we see that the density of
quasiparticles is bounded approximately as $n\sim e^{-\beta\Delta}$. The
apparent saturation of MCP around $\beta\Delta=4\pi$ is due to small ringing
errors, which can be removed by increasing the ratio $T/\beta$ — we refer the
reader to Appendix A for more details. In practice, dominant heating processes
arise due to noise (see Sec. IV) and it pays to keep the ratio $T/\beta\gtrsim
O(1)$.

From the above results, we conclude that the MCP is able to successfully
prepare the ground state of the TFIM, with fidelity approaching 1 in the limit
of large $\beta$. However, the associated time scale diverges with the system
size, naturally a weakness of the edge cooling setup. The rate equation
accurately captures the quasiparticle dynamics for the case of edge cooling.
Next, we will consider the bulk cooling setup with a finite density of
auxiliaries.

###  III.3 Finite density of auxiliaries

Figure 3: Performance of the bulk cooling protocol, with one auxiliary per
system qubit and $N_{S}=20$. The parameter choice is described in the main
text. (a) Quasiparticle density $n(t)$ vs. rescaled timescale
$t=\nu\theta^{2}$ in the PM phase. The MPS results (solid lines) are compared
to the kinetic theory predictions (dashed lines), for both SCP and MCP.
Parameter $\theta=0.1$. (b) Same as in (a) but for the AFM phase. We show MPS
results for two coupling strengths, $\theta=0.1$ and $0.05$, with the larger
showing corrections to the rate equation. The AFM cools more slowly than the
PM, as expected due to the topological nature of the quasiparticles. (c) Log-
fidelity per qubit of the steady state with quasiparticle vacuum, calculated
from the kinetic equation, as a function of parameter $\beta$, related to the
number of unitary layers in one cooling cycle, via $T=5\beta$.

While providing an instructive example, the edge cooling setup is not
practical for large systems. We argued above that the cooling timescale
diverges with the system size in that case, which holds generally for models
with local boundary dissipation [59]. Furthermore, the cooling rate remains
constant while the heating rate due to noise scales with system size (see the
following Section). Therefore, we now consider the setup in Fig. 1d, with a
finite density of auxiliaries. For simplicity, we will assume that there is
one auxiliary per system qubit.

In this setup, the dynamics is no longer integrable, and the system’s density
matrix develops non-Gaussian correlations due to the fact that the matrix
elements for bulk $\sigma_{j}^{\pm}$ operators carry Jordan-Wigner strings in
the fermionic formulation (see Eq. (32)). This allows for transitions between
many-body eigenstates where the occupation numbers of several quasiparticle
levels change at once.

To make progress, we therefore assume that the processes involving a small
number quasiparticles are dominant. For cooling processes, this assumption is
justified in the low-density limit $n_{k}\ll 1$, which describes the system’s
dynamics close to the ground state 444Heating processes, in particular due to
noise, may in principle involve a large number of particles; however, we will
find that few-quasiparticle processes are dominant, and considering them gives
an accurate approximation of heating rates.. Restricting to processes
involving at most two quasiparticles, we derive the following kinetic
equation, containing cooling and heating terms, as well as terms that describe
scattering between quasiparticles:

$\displaystyle\begin{split}\delta
n_{k}&=-n_{k}W_{k}^{-}+(1-n_{k})W_{k}^{+}+\sum_{q}\Big{[}-n_{k}W_{k,q}^{-}n_{q}+(1-n_{k})W_{k,q}^{+}(1-n_{q})-n_{k}V_{k,q}^{-}(1-n_{q})+(1-n_{k})V_{k,q}^{+}n_{q}\Big{]}.\end{split}$
(41)

We compute the rates in Eq. (41) assuming the system is near the vacuum state
$\Omega$, which yields:

$W_{k}^{\mp}=\theta^{2}\sum_{j=1}^{N_{S}}|F_{h,T}(\pm\epsilon_{k})|^{2}|\sigma^{\pm}_{j;0;k}|^{2},$
(42)
$W_{k,q}^{\mp}=\theta^{2}\sum_{j=1}^{N_{S}}|F_{h,T}(\pm(\epsilon_{k}+\epsilon_{q}))|^{2}|\sigma^{\pm}_{j;0;k,q}|^{2},$
(43)
$V_{k,q}^{\mp}=\theta^{2}\sum_{j=1}^{N_{S}}|F_{h,T}(\pm(\epsilon_{k}-\epsilon_{q}))|^{2}|\sigma^{\pm}_{j;q;k}|^{2},$
(44)

where we introduced the shorthand matrix-elements

$\sigma^{\pm}_{j;k_{1},\ldots,k_{m};k_{1},\ldots,k_{n}}=\bra{k_{1},\ldots,k_{m}}\hat{\sigma}^{\pm}_{j}\ket{k_{1},\ldots,k_{n}}.$
(45)

Above, the terms $W^{\mp}_{k}$ correspond to cooling/heating processes
involving one quasiparticle, while the terms $W^{\mp}_{k,q}$ describe
removing/adding a pair of quasiparticles, and $V^{\mp}_{k,q}$ describe
quasiparticle scattering. Note that while the scattering does not directly
change the total quasiparticle occupation number, it will tend to move
quasiparticles toward lower quasienergies due to the presence of the filter
function $F_{h,T}$. We expect that these rates will be weakly affected by a
small finite background density of quasiparticles. The matrix elements can be
evaluated using Wick’s theorem (see Appendix D and a review [61]).

Figure 4: Quasiparticle occupation numbers $n^{\infty}_{k}$ in the steady
state, as a function of quasienergy $\epsilon_{k}$. We compare MPS simulations
(scatter points) and kinetic equation results (red line), for a system with
$N_{S}=N_{A}=20$ qubits in (a) PM and (b) AFM phases. The reset time is
$T=12=2\beta$, and the coupling $\theta=0.05$. The MPS simulations use a
maximum bond dimension of $d=100$. Inset: deviation $|\delta n_{k}|$ between
rate equation and MPS values shown in main figure.

Crucially, the matrix elements entering in the above rates obey different
selection rules in the PM and AFM phases. In the PM phase, there is a unique
symmetric ground state. Operators $\hat{\sigma}_{j}^{\pm}$ anticommute with
the parity operator $\prod_{j}\hat{Z}_{j}$, so in the fermionic language they
must change the fermion parity. Thus, two particle processes are forbidden in
the PM phase, $W^{\mp}_{k,q}=V^{\mp}_{k,q}=0$. In contrast, in the AFM phase,
there are two vacuum states, $|\Omega_{\pm}\rangle$, which have parity $\pm
1$. The cooling process can induce transitions between them by changing the
occupation of the Majorana edge fermion. $W_{k,q}^{\mp}$ and $V_{k,q}^{\mp}$
matrix elements are non-vanishing in this case. The bra and ket states should
be understood as quasiparticle states on top of different vacua. Note that the
Majorana edge mode has a vanishing energy in the limit of large system size;
hence, this mode cannot by cooled by our protocol and will remain effectively
at an infinite temperature, corresponding to an equal weight mixture of the
two vacuum states.

In the spin language, these selection rules follow from intuitive arguments
given in Ref. [43]. The quasiparticles of the PM phase are magnons, local
excitations that can be removed via local spin flips (caused by auxiliary
qubit interactions); in contrast, in the AFM phase, quasiparticles are domain-
wall-like and can only be removed/created locally in pairs. In Appendix D, we
study the distinct behavior of matrix elements in the two phases in detail. In
particular, the two-particle matrix elements in the AFM phase are homogeneous
in the bulk, but single-particle matrix elements decay exponentially with the
distance of the auxiliary qubit from the boundary (since single domain walls
can be removed at the edges of the chain), with a decay length set by the
correlation length of the model. We will see below that this difference in
quasiparticle nature leads to different efficiencies of the cooling in the two
phases.

We now compare the predictions of the above kinetic theory to matrix product
state (MPS) simulations of the system’s density operator dynamics (more
precisely, density matrix product operator simulations [62, 63]). The dynamics
are performed using the standard time-evolving block decimation (TEBD)
algorithm [64, 65], implemented via the ITensor library [66]. Implementation
details and numerical verifications of asymptotic convergence are presented in
Appendix E.

As illustrated in Fig. 3a,b, the theory continues to be accurate for cooling
with a finite density of auxiliaries, for both the SCP and MCP protocols. We
fix system size $N_{S}=20$ with one auxiliary per site. We use cooling
parameters $T=5$, $h=0.3$ in the SCP and $T=12=2\beta$ in the MCP. The
coupling is chosen to be $\theta=0.1$ in the PM phase. In the AFM phase, this
relatively large coupling leads to visible deviations between the MPS results
and the kinetic equation prediction at late times; hence we include simulation
results for a weaker coupling, $\theta=0.05$, which leads to a good agreement.
We use bond dimension $d=100$ for both the SCP and MCP ($\theta=0.05$) data,
and $d=150$ for the MCP ($\theta=0.1$) data. These choices are converged in
bond dimension (see Appendix E). Our choice of larger $\theta$ and smaller MCP
reset time (compared to the edge cooling case) helps to avoid high
entanglement at intermediate times, which affects the accuracy of the MPS
simulations. The MCP outperforms the SCP, and in Fig. 3c we provide further
confirmation that the MCP can reach much higher fidelities, plotting the
increase in the log-fidelity per qubit as the parameter $T=5\beta$ is
increased (we choose a large value of $T/\beta$ to remove the ringing
artifacts for the range of fidelities shown, see Appendix A).

The difference in cooling distinct quasiparticle types in the PM and AFM
phases is also evident in the timescales shown in Fig. 3a,b: the PM reaches a
quasiparticle density below $10^{-3}$ after a small number of cooling cycles
$O(\theta^{-2})$, with a clear exponential decay to the steady state value. In
contrast, in the AFM phase, cooling slows down once quasiparticle density is
below $n(t)\lesssim 10^{-1}$. This separation of timescales in the AFM is due
to the fact that initial cooling proceeds quickly through pair annihilation
processes, but this mechanism is suppressed $\propto n^{2}$ in the low density
regime. Late-time cooling is therefore dominated by single-particle processes,
which may only occur around the system edge.

As a final comparison between the kinetic equation and the MPS numerics, we
look at the momentum-resolved quasiparticle occupations in the steady state,
in Fig. 4, focusing on the MCP. To reach the steady state in our MPS
simulations, we perform the simulations starting from the exact ground state.
As we expect the (unique) steady state to be close to the true ground state,
this allows for faster real-time simulations compared to starting from a high
energy-density state. The figure illustrates that occupation numbers of all
levels are below $10^{-2}$, and the remaining population is mostly due to
quasiparticles at the band edge. The finite weight on low energy modes stems
from the relatively small value of the parameter $\beta\approx
1.2\pi\Delta^{-1}$ chosen to aid the convergence of the MPS simulations.

To sum up, the analysis above demonstrates efficiency of MCP with a finite
density of auxiliaries for removing quasiparticles after $O(\theta^{-2})$
cooling cycles, in the PM phase. In agreement with the arguments of Ref. [43],
we also find that topological, domain-wall type quasiparticles are more
challenging to remove. Nevertheless, in an idealized system without noise, our
cooling protocol allows preparation of the ground state with fidelity
arbitrarily close to 1.

##  IV Noise effects

Having established that the quasiparticle cooling algorithm can be used to
efficiently prepare correlated ground states, we now turn to the effects of
unwanted decoherence. Each experimental platform comes with its own sources of
noise, and our goal will be to understand the limitations it imposes on the
cooling efficiency. We will focus on the single-qubit decay and dephasing
processes, relevant in systems of superconducting qubits. The kinetic equation
can be generalized to include transitions mediated by weak noise. Using this
modified kinetic theory, we analyze the dependence of the steady state
quasiparticle population on the noise strength.

###  IV.1 Rate equation in the presence of noise

Considering continuous time evolution, single-qubit decoherence processes can
be conveniently described within Lindblad formalism [50]. The evolution of the
system’s density matrix $\rho_{S}$ is given by (for the moment, we ignore
presence of auxiliary qubits):

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=-i[\hat{H},\hat{\rho}_{S}]+\delta{\mathcal{D}}(\hat{\rho}_{S};\tilde{\gamma}),$
(46)
$\displaystyle\delta{\mathcal{D}}(\hat{\rho}_{S};\tilde{\gamma})=\sum_{j,\mu}\tilde{\gamma}_{j\mu}\left(\hat{L}_{j\mu}\hat{\rho}_{S}\hat{L}_{j\mu}^{\dagger}-\frac{1}{2}\\{\hat{L}_{j\mu}^{\dagger}\hat{L}_{j\mu},\hat{\rho}_{S}\\}\right),$

where index $\mu$ runs over possible decoherence processes,
$\tilde{\gamma}_{j\mu}$ are their rates, and $\hat{L}_{j\mu}$ are the
corresponding jump operators.

In the context of Floquet evolution, we will adopt the following simplified
model: we assume that decoherence can be described by a quantum channel
$\mathcal{D}$ acting on the system after one period of Floquet evolution. The
quantum channel will be constructed by considering only the dissipative part
$\delta\tilde{\mathcal{D}}$ of the Lindblad generator in Eq. (46), acting over
the Floquet cycle duration $\Delta t$. Experimentally, quantum processors
operate in the regime where $\tilde{\gamma}_{j\mu}\Delta t\ll 1$; then, this
dissipative channel can be written as (to leading order in $\tilde{\gamma}$
coefficients):

$\mathcal{D}(\hat{\rho}_{S})\approx\hat{\rho}_{S}+\delta\mathcal{D}(\hat{\rho}_{S};\gamma),$
(47)

where we introduced dimensionless decoherence rates:

$\gamma_{j\mu}=\tilde{\gamma}_{j\mu}\Delta t,\;\;\gamma_{j\mu}\ll 1.$ (48)

While this approximate model is only expected to be quantitatively accurate in
the Trotter limit $g,J\ll 1$, we expect it to correctly capture the
qualitative effects of noise on cooling even away from that limit. The main
motivation for introducing this simplified model is that it decouples the
unitary evolution of the cooling protocol from the decoherence effects,
leading to a transparent and compact form of the noisy kinetic equation.

Within this model, the evolution of the system-bath density matrix over one
cooling period is given by:

$\Phi^{\gamma}(\rho_{S})=\text{Tr}_{B}\
\mathcal{D}\mathcal{M}_{T}...\mathcal{D}\mathcal{M}_{1}(\hat{\rho}_{S}\otimes\hat{\rho}_{B}^{0}),$
(49)

where

$\mathcal{M}_{\tau}(\cdot)=\hat{U}_{\theta,\tau}\hat{U}_{B}\hat{U}_{S}(\cdot)\hat{U}_{S}^{\dagger}\hat{U}^{\dagger}_{B}\hat{U}_{\theta,\tau}^{\dagger}.$

We follow the same steps that led to Eq. (23) above, now keeping first-order
terms in $\theta^{2}$ and $\gamma_{j\mu}$. The processes from cooling and
noise contribute independently, since we discard next-order terms
$O(\gamma\theta^{2})$. The change in the quasiparticle occupations due to
noise takes the standard Lindblad form

$\delta\hat{n}_{k}=T\sum_{j,\mu}\gamma_{j\mu}\left(\hat{L}^{\dagger}_{j\mu}\hat{n}_{k}\hat{L}_{j\mu}-\frac{1}{2}\\{\hat{L}_{j\mu}^{\dagger}\hat{L}_{j\mu},\hat{n}_{k}\\}\right),$
(50)

and following Eqs. (16-23) we obtain the noisy rate equation, which includes
cooling transitions $R_{\vec{\alpha}\vec{\beta}}$ and and an additional
contribution $R_{\vec{\alpha}\vec{\beta}}^{\gamma}$ to transition rates from
the noise terms:

$\displaystyle\delta n_{k}$
$\displaystyle=\sum_{\vec{\alpha},\vec{\beta}}\rho_{\vec{\alpha}\vec{\alpha}}(\beta_{k}-\alpha_{k})[R_{\vec{\alpha}\vec{\beta}}+R_{\vec{\alpha}\vec{\beta}}^{\gamma}],$
(51) $\displaystyle R_{\vec{\alpha}\vec{\beta}}$
$\displaystyle=\theta^{2}\sum_{j=1}^{N_{A}}|F_{h,T}(\Delta(\vec{\alpha},\vec{\beta}))|^{2}|\bra{\vec{\beta}}\hat{\sigma}_{Q_{j}}^{+}\ket{\vec{\alpha}}|^{2},$
(52) $\displaystyle R_{\vec{\alpha}\vec{\beta}}^{\gamma}$
$\displaystyle=\sum_{j=1}^{N_{S}}\sum_{\mu}\gamma_{j\mu}T|\langle\vec{\beta}|\hat{L}_{j\mu}|\vec{\alpha}\rangle|^{2}.$
(53)

Note that the noise is essentially ‘infinite temperature’ and has no term
analogous to the spectral filtering function $F_{h,T}$.

###  IV.2 Application: cooling the transverse-field Ising model

Figure 5: Performance of the cooling algorithm under decoherence. (a)
Quasiparticle density in the steady state, as a function of noise ratio
$\gamma/\theta^{2}$, in PM (red) and AFM (blue) phases, for different system
sizes $N_{S}$ from 10 to 80. We use MCP with fixed $T=2\beta=6\pi/\Delta$. At
weak noise, PM quasiparticle density scales linearly with $\gamma$, while AFM
displays a crossover from $n\propto\gamma$ to $n\propto\sqrt{\gamma}$
behaviour as noise strength is reduced. (b) Log-fidelity between steady state
and quasiparticle vacuum, for MCP (solid lines) and SCP (dashed lines). We fix
$T=2\beta$ and optimize over parameter $\beta$ for each protocol.

Next, we apply the noisy kinetic equation to cooling of a transverse-field
Ising model in the presence of decay to the qubit ground state, and dephasing
noise. This corresponds to $\mu=d,\phi$, with jump operators

$\hat{L}_{jd}=\hat{\sigma}_{j}^{+},\hskip
14.22636pt\hat{L}_{j\phi}=\hat{Z}_{j}.$ (54)

We assume homogeneous decoherence rates for simplicity,
$\gamma_{j\phi}=\gamma_{\phi}$, $\gamma_{jd}=\gamma_{d}$, and rescale the
dephasing rate, $\gamma_{\phi}\to\gamma_{\phi}/2$, following standard
conventions in the literature on open quantum system.

Similar to our analysis in Sec. III, we restrict to processes involving at
most two quasiparticle levels. (Processes involving a higher number of
quasiparticles can be systematically accounted for within our formalism.) This
approximation yields the same kinetic equation derived in Sec. III, Eq. (41),
but with the rates modified by noise, $W\to\tilde{W}$, $V\to\tilde{V}$:

$\displaystyle\tilde{W}_{k}^{\mp}={W}_{k}^{\mp}+\sum_{j=1}^{N_{S}}\gamma_{d}T|\sigma^{\pm}_{j;0;k}|^{2}+\frac{\gamma_{\phi}T}{2}|Z_{j;0;k}|^{2},$
(55)
$\displaystyle\tilde{W}_{k,q}^{\mp}=W_{k,q}^{\mp}+\sum_{j=1}^{N_{S}}\gamma_{d}T|\sigma^{\pm}_{j;0;k,q}|^{2}+\frac{\gamma_{\phi}T}{2}|Z_{j;0;k,q}|^{2},$
(56)
$\displaystyle\tilde{V}_{k,q}^{\mp}=V_{k,q}^{\mp}+\sum_{j=1}^{N_{S}}\gamma_{d}T|\sigma^{\pm}_{j;q;k}|^{2}+\frac{\gamma_{\phi}T}{2}|Z_{j;q;k}|^{2}.$
(57)

The matrix elements satisfy the selection rules discussed in Sec. III. We now
summarise the main behaviour of the noisy rate equation in the two phases,
with a more detailed discussion of matrix elements for $\hat{Z}_{j}$ and
$\hat{\sigma}^{\pm}_{j}$ operators given in Appendix D. In the PM phase,
single particle loss and gain processes occur throughout the chain (cooling
and decay contributions); quasiparticle pair creation/annihilation and number-
conserving scattering events occur due to dephasing. In the AFM phase, on
account of the domain-wall nature of excitations, only two-particle scattering
and pair creation/annihilation processes can occur in the bulk, and they are
induced by cooling, decay and dephasing. Single particle transitions require
an interaction with the Majorana edge mode, and are localized to the chain
edge.

We have compared the predictions of the noisy kinetic equation and MPS
simulations in the presence of noise, finding again a good agreement in the
parameter regimes accessible to MPS (moderate value of noise strength, or
moderate times at a very weak noise). We refer the reader to Appendix F for a
comparison of the steady state occupations, also for different coupling
strengths. The rate equation, however, allows us to study system sizes and
parameter regimes far beyond the reach of matrix product state simulations
555The steady state can be efficiently found by root-finding algorithms, after
the cooling and noise-induced rates have been calculated via Wick’s theorem.
In practice, the system reaches its steady state quickly and we evolve the
rate equation in time until the quasiparticle occupations have effectively
converged..

In Fig. 5a, we illustrate the quasiparticle density in the steady state,
obtained from the kinetic equation, as a function of the noise strength. For
simplicity, we choose $\gamma_{d}=\gamma_{\phi}=\gamma$ and fix auxiliary
density $N_{A}=N_{S}$. We focus on the MCP with $T=2\beta=6\pi/\Delta$. We
observe a qualitatively different behavior in the two phases, in line with
above results on noiseless cooling. In the PM phase, the steady state
quasiparticle density decays as a power-law function of noise,
$n^{\infty}_{\rm PM}\propto\gamma^{a}$, $a\approx 1$, with little dependence
on the system size. In the AFM phase, in contrast, the density displays a
power-law scaling which depends on system size, $n^{\infty}_{\rm
AFM}\propto\gamma^{b(N_{S})}$. For small systems and weak noise, we observe
similar scaling to the PM, with $b\approx 1$; as either the system size or
noise strength is increased, we find a crossover to a different power
$b\approx 1/2$.

The utility of the kinetic theory is demonstrated with the following
simplified versions of the kinetic equations, which capture the observed
features in Fig. 5a: In the PM phase, where single-quasiparticle processes can
occur throughout the system (bulk and edge), the evolution of the occupation
number $n_{k}(t)$ is given by

$\frac{dn_{k}(t)}{dt}\approx-C_{1}\theta^{2}n_{k}+C_{\gamma}\gamma T,$ (58)

where $C_{1},C_{\gamma}$ are constants of order one. We assumed that MCP is in
the regime where heating processes due to coupling to auxiliaries are strongly
suppressed, and neglected the dependence of $C_{1},C_{\gamma}$ on $k$. For the
noise we neglect terms $O(\gamma Tn)$, since we are interested in the regime
of weak noise and low density. Then, the expected quasiparticle density
$n^{\infty}$ is given by

$n^{\infty}_{\rm PM}\approx\frac{C_{\gamma}\gamma T}{C_{1}\theta^{2}},$ (59)

where we averaged over $k$, which gives the linear dependence on $\gamma
T/\theta^{2}$ in Fig. 5a, for weak noise.

In the AFM phase, as discussed above, the single-quasiparticle cooling
processes can only occur at the edge, and their rates are suppressed by

$\theta^{2}\Gamma_{k}=\theta^{2}C_{e}|u_{1,k}|^{2}\approx\theta^{2}C_{e}\frac{k^{2}}{N_{S}}$
(60)

(strictly speaking, this expression holds for small quasimomenta, $k\ll 1$,
but we will see below that the steady state quasiparticle density is
determined mostly by this quasimomentum range). The processes removing
quasiparticle pairs can be modeled by a term $-C_{2}\theta^{2}n_{k}n$,
neglecting the $k,q$ dependence of the rates $W_{k,q}^{\mp}$. The dominant
heating mechanism is from pair creation in the bulk, $C_{\gamma}\gamma
T(1-n_{k})(1-n)$. Then, we arrive at

$\displaystyle\frac{dn_{k}(t)}{dt}\approx-\theta^{2}\Gamma_{k}n_{k}-C_{2}\theta^{2}n_{k}n$
$\displaystyle+C_{\gamma}^{\prime}\gamma T(1-n_{k})(1-n).$ (61)

The steady state occupation numbers are given by

$n_{k}^{\infty}\approx\frac{C_{\gamma}^{\prime}\gamma
T}{\theta^{2}\Gamma_{k}+C_{2}\theta^{2}n^{\infty}_{\rm
AFM}+C_{\gamma}^{\prime}\gamma T},$ (62)

supplemented by a self-consistency condition

$n^{\infty}_{\rm AFM}=N^{-1}_{S}\sum_{m=1}^{N_{S}}n^{\infty}_{k_{m}},$ (63)

where we again neglected terms $O(\gamma Tn)$ for the regime of interest. The
solution $n^{\infty}_{\rm AFM}(\gamma)$ exhibits two regimes, depending on the
noise strength. At stronger noise values, where $n^{\infty}N_{S}\gg 1$ (such
that the total number of quasiparticles in the steady state of the system is
much greater than 1), two-particle cooling processes dominate. Neglecting the
sub-dominant single-particle cooling contribution $\Gamma_{k}$ yields

$n^{\infty}_{\rm AFM}\approx\left(\frac{C_{\gamma}^{\prime}\gamma
T}{C_{2}\theta^{2}}\right)^{1/2}.$ (64)

In the opposite limit, $n^{\infty}N_{S}\ll 1$, the single-particle cooling
processes determine the steady state and we can neglect the $C_{2}$ term in
Eq. (62). Then,

$\displaystyle n^{\infty}_{\rm
AFM}\approx\sum_{m=1}^{N_{S}}\frac{C_{\gamma}^{\prime}\gamma
T}{C_{e}\theta^{2}m^{2}N_{S}^{-2}+C_{\gamma}^{\prime}\gamma TN_{S}}\approx
C\frac{\gamma TN_{S}^{2}}{\theta^{2}},$ (65)

with $C=\pi^{2}C_{\gamma}^{\prime}/6C_{e}$. This expression holds in the weak-
noise limit, $n^{\infty}N_{S}\sim{\gamma TN_{S}^{3}}/{\theta^{2}}\ll 1$, and
we made use of the approximate quantization condition
$k_{m}\approx\pi(1-m)/N_{S}$, which holds for $k_{m}\ll\pi$. These approximate
rate equations, Eqs. (64,65), capture the crossover behaviour displayed in
Fig. 5a.

In Sec. III we argued that the fidelity of the steady state with the ground
state approaches 1 as the MCP parameter $T\propto\beta$ is made large. In the
presence of noise, there will be an optimal $T_{*}$ and $\beta_{*}$ that
maximises the fidelity, since the noise strength also scales with $T$. We
compare the log-fidelity per qubit for the MCP and the SCP vs. the noise ratio
$\gamma/\theta^{2}$, in Fig. 5b. For each point, we optimize the parameter
$\beta$, fixing $T=2\beta$ for simplicity (for the SCP this simply corresponds
to optimizing $T$). We observe that, while at moderate noise values
$\gamma/\theta^{2}\gtrsim 10^{-2}$, the SCP marginally outperforms MCP, for
weak noise the MCP reaches lower fidelities in both phases, approaching
fidelity 1 as the noiseless limit is approached. We conclude that the MCP is
able to prepare target states with high fidelity even in the presence of weak
noise, but noise is more destructive in systems with topological excitations
such as the 1d AFM example.

##  V Cooling non-integrable systems

In this Section, we investigate the efficiency of the MCP for preparing ground
states of non-integrable models. While these models are generally
thermalizing, they typically have long-lived quasiparticles at low energy
density. Therefore, the theory of dissipative cooling described in the
previous sections for integrable systems, is expected to hold, once the system
has been cooled to a state with sufficiently low energy density.

For a numerical test, we focus on the example of the TFIM with an additional
longitudinal field,

$\displaystyle\hat{U}^{\prime}_{S}=e^{-i\hat{G}_{0}}e^{-i\hat{G}_{1}},$ (66)
$\displaystyle\hat{G}_{0}=\frac{\pi
J}{2}\sum_{i=1}^{N_{S}-1}\hat{X}_{i}\hat{X}_{i+1},\quad\hat{G}_{1}=-\frac{\pi}{2}\sum_{i=1}^{N_{S}}g\hat{Z}_{i}+g_{X}\hat{X}_{i}.$

We fix $J=0.1,g=g_{X}=0.15$, for which the model is in a paramagnetic phase
[48]. As discussed above, since quasiparticles in the PM phase are simple
magnons, we expect this to be the fastest regime to cool, and simulation
accuracy is higher since the system reaches the low-density regime much
faster.

Our parameter choice is in the Floquet regime, and not in the strict Trotter
limit ($J,g_{x},g_{z}\ll O(1)$). Thus, to characterize the cooling efficiency,
we consider effective Hamiltonians obtained in different orders of the
Baker–Campbell–Hausdorff (BCH) expansion of Eq. (66),
$\hat{U}_{S}^{\prime}\approx e^{-i\hat{H}_{l}}$, $l=0,1,\ldots$ (we will limit
our consideration to $l=0,1$.) The first two orders of BCH formula give

$\hat{H}_{0}=\hat{G}_{0}+\hat{G}_{1},\quad\hat{H}_{1}=\hat{G}_{0}+\hat{G}_{1}-\frac{i}{2}[\hat{G}_{0},\hat{G}_{1}].$
(67)

The cooling efficiency can then be characterized by measuring the energy and
the fidelity of the system’s density matrix with respect to the ground states
of the above Hamiltonians:

$E_{l}=\text{tr}(\hat{H}_{l}\hat{\rho}),\qquad\mathcal{F}_{l}=\langle\Omega_{l}|\hat{\rho}|\Omega_{l}\rangle,$
(68)

where $l\in\\{0,1\\}$ denotes the BCH expansion order, and $\Omega$ is the
corresponding ground state. The system is in the prethermal regime, due to the
fact that $J,g,g_{X}\sim O(10^{-1})$, and is expected to cool towards the
ground state of an effective Hamiltonian $\hat{H}_{l_{*}}$, with an optimal
truncation order $l_{*}\geq 1$ of the BCH series. Thus, we expect the ground
state of $\hat{H}_{1}$ to provide a better approximation than that of
$\hat{H}_{0}$.

Figure 6: Cooling a non-integrable spin chain. (a) Evolution of energy,
defined with respect to the ground state of approximate effective Hamiltonian
$H_{0}$ (dashed lines) and $H_{1}$ (full lines), as a function of the number
of resets $\nu$, for different system sizes. Corresponding energy values for
the first excited states of $H_{1}$ are indicated by arrows, showing
convergence of energy to a value well below the excitation gap. (b) Log-
fidelity per qubit vs. number of resets. Model parameters are
$h_{x}=h_{z}=0.15,J=0.1$ (PM phase), reset time $T=20=2\beta$. The bond
dimension of simulations is $d=300$.

For the cooling protocol we focus on the MCP with finite density
$N_{A}=N_{S}$. In contrast to the previous sections, where our aim was to
develop a theory of cooling, here we generalise the protocol to cool more
rapidly by using a large initial value of $\theta$, later approaching the weak
coupling limit $\theta\ll 1$ where the previous analysis is expected to hold.
To this end, we initialize the coupling at $\theta(0)\approx 0.4$ and
progressively reduce to a final value $\theta_{\infty}\approx 0.083$. The
value of the coupling is reduced as $\theta(t+1)=0.8\theta(t)$ whenever the
energy of the system is reduced slower than a chosen rate,
$0<E_{1}(t)-E_{1}(t-1)<0.01$. We found this choice to be efficient for the
finite time simulations. Note that the function $f_{\tau}$ is still modulated
within a cycle, see Eq. (29). We choose a reset time $T=20=2\beta$, for which
$\Delta\beta\approx 2.6$, where $\Delta$ is the (many-body) energy gap of
$\hat{H}_{1}$.

The tensor-network simulations are performed similarly to the integrable case,
see Appendix E for details. The ground and first excited eigenstates of
$\hat{H}_{0},\ \hat{H}_{1}$, and their corresponding energies, used in (68),
are evaluated with a standard DMRG algorithm; the eigenstates are truncated to
a bond dimension $d_{0}\approx\sqrt{d}$, where $d$ is the maximal bond
dimension of the MPS performing the cooling protocol.

In Fig. 6a, we show how the energies defined with respect to $\hat{H}_{0}$ and
$\hat{H}_{1}$ evolve over the course of the cooling protocol, for three
different system sizes $N_{S}=10,\ 20,\ 40$. We observe that the system energy
decays quickly toward the ground state energy of $\hat{H}_{1}$. (Note that due
to the sweeping of $\theta$ the decay rate varies in time, as opposed to the
clean exponential decay in the integrable case.) The expectation value of the
zero-order Hamiltonian $\hat{H}_{0}$ saturates at a higher value than that of
$\hat{H}_{1}$, indicating that the steady-state is well-approximated by the
ground state of $\hat{H}_{1}$, as expected. Cooling efficiency appears to
slightly increase with larger system sizes at later times, but larger system
studies are required to provide a conclusive answer. In Fig. 6b, we observe
similar behaviour in the decay of the log-fidelity per qubit, further
corroborating the ability of MCP to cool the system to the ground state of an
effective prethermal Hamiltonian.

These results show that the MCP remains an effective tool for cooling non-
integrable systems, which exhibit long-lived quasiparticles at a sufficiently
low energy density. We expect that the performance of cooling protocols can be
further improved, e.g. by adapting methods employed in variational quantum
algorithms [68].

##  VI Discussion and outlook

In summary, in this work we investigated cooling protocols for preparing
ground states of many-body systems on quantum processors. Our approach is
reminiscent of cooling strategies in systems of ultracold atoms, including
sympathethic cooling [69], but our protocol primarily aims to remove
quasiparticles – emergent excitations above a many-body ground state. We
developed a phenomenological kinetic theory framework, which accounts for
competing cooling and heating processes induced by auxiliaries and by noise,
and compared the efficiency of different cooling protocols.

Our results are directly applicable to studying correlated states with current
digital quantum processors, as well as in analog processors with the fast
reset capability: indeed, SCP was recently implemented on the Google Quantum
Processor [38]. While reaching the ground state in that experiment was
challenging due to the level of noise $\gamma\approx 0.01$, our work shows
that, with e.g. $\theta=0.1$, noise values of $\gamma\approx 10^{-3}$ will be
sufficient for preparing steady states with quasiparticle density below
$n\approx 10^{-2}$ in the paramagnetic phase (see Fig. 5b). In systems of
$N_{S}=100$ qubits, this translates to approximately one quasiparticle in the
sample. Cooling the antiferromagnetic phase is more challenging; nevertheless,
noise level $\gamma\approx 10^{-6}$ allows to reach comparable densities. We
expect that such low-energy states can be further efficiently purified to
extract properties of the ground state and low-energy dynamics [38].

The quasiparticle cooling algorithm is a promising alternative to existing
state-preparation methods, including variational unitary circuits [19, 20,
21], thanks to a degree of robustness with respect to noise it offers, as well
as its generality. First, dissipation removes excitations created by noise,
and the steady state is independent of the initial state and the history of
errors. Second, our study of non-integrable TFIM suggests that the cooling
approach may be suitable for preparing a broad class of correlated states in
non-integrable systems. We note that the cooling algorithm requires the
ability to engineer a local coupling operator which has a non-zero overlap
with quasiparticle operators. Developing strategies to identify efficient
auxiliary-system coupling for models of quantum magnetism, in particular in
2d, is an important direction for future research.

Despite the generality of this approach, we expect certain states will yet
prove hard to prepare. One reason is that preparing ground states with high
fidelity in general would allow one to solve QMA-complete problems [70], which
would contradict a widely believed conjecture that QMA-complete problems
cannot be solved efficiently on quantum computers. Still, we do not expect the
extreme hard cases to arise in the absence of spin-glass order. In this
context, it is therefore important to understand which phases of matter can be
prepared by the quasiparticle cooling and other dissipative algorithms.

As pointed out in Ref. [43], and further highlighted by our analysis, systems
with topological, non-local quasiparticles are more challenging to cool. On
the fundamental level, this stems from lower bounds [71, 72] on the
preparation times of topologically ordered states (e.g., toric code).
Exploring strategies of optimal cooling of such systems, for example by
creating quasiparticle traps or making use of dynamical gauge fields [73], is
a promising direction of future research. Furthermore, it would be interesting
to apply our analysis to quantum-critical points with gapless quasiparticles,
and to study the possibility of cooling strongly interacting systems without
long-lived quasiparticles. Generally, we may expect [74] that algorithms
similar to ours will be able to efficiently prepare low-energy states of
_thermalizing systems_ (which rules out examples such as spin glasses, where
cooling procedure will drive the system to one of the many metastable minima).

Our work suggests several other avenues for future research. First, we note
that our protocol may be modified to prepare approximate finite-temperature
Gibbs states [27, 74, 75, 34], by appropriately choosing the filter function
and the operator entering coupling to auxiliary, such that the corresponding
rates in Eqs. (42-44) approximately satisfy the detailed balance condition. In
fact, this is achieved by our MCP, Eq. (29), see Appendix B, giving a
corresponding effective temperature $T_{\rm eff}=1/2\beta$. Weak noise
processes in quantum processors generally violate detailed balance, and will
make the steady state slightly non-thermal. In contrast, intrinsic
thermalization mechanisms in interacting systems are expected to favor the
thermal state. This competition, and the resulting steady state may be studied
using digital [38] and analog quantum processors.

Finally, our protocol can be modified to engineer non-equilibrium
quasiparticle distribution in the steady state. The properties of the steady
state, and evolution towards it can give useful insights into interactions
between quasipaticles, their dynamics, and lifetimes – properties that play a
central role in our understanding of correlated materials. We note that the
accuracy of the kinetic theory we used to model the cooling protocol suggests
that the dominant thermalization mechanisms are well approximated by few-
particle processes. Applications of suitably modified protocols in quantum
simulator experiments would unlock new avenues to study underlying dynamics
and thermalization mechanisms.

###### Acknowledgements.

We thank Zala Lenarčič, Iris Ulčakar and Guifre Vidal for insightful
discussions, and Rolando Somma and Guifre Vidal for useful comments on the
manuscript. This work was partially supported by the European Research Council
via the grant agreements TANQ 864597 (A.M., D.A.).

##  Appendix A Details of modulated coupling pulse

In this Appendix, we provide additional details regarding the modulated
coupling protocol (MCP) defined in Eq. (29). As explained in Section II, this
pulse is advantageous in that the corresponding filter function,
$F_{h,T}(\epsilon)$, approximates a broadened step function, with broadening
controlled by the parameter $\beta$. This can be seen as follows: Recalling
the definition of the filter function (Eq. 24), we have for the MCP

$F_{h,T}(\epsilon)=\frac{\pi}{A\beta}\sum_{\tau=1}^{T}\frac{\sin[\pi(\tau-\tau_{0})/2]}{\sinh[\pi(\tau-\tau_{0})/\beta]}e^{i\tau(\epsilon-\pi
h)}.$ (69)

For the MCP we choose the auxiliary field $h=1/2$. We consider the limits
$T,\beta\to\infty$, with the ratio $\alpha=\beta/T\ll 1$ held constant. Then,
the summation (69) can be replaced with the integral

$\frac{\pi}{A}\int_{-1/2\alpha}^{1/2\alpha}dz\ \frac{e^{i\epsilon\beta
z}-e^{i(\epsilon-\pi)\beta z}}{2\sinh(\pi z)},$ (70)

up to an unimportant phase factor. For $\alpha\ll 1$ the integral limits can
be extended to infinity and (70) evaluated via residues, giving ($\epsilon>0$)

$F_{h,T}(\epsilon)=\pi\bigg{[}\frac{\tanh(\epsilon\beta/2)-\tanh((\epsilon-\pi)\beta/2)}{\tanh(\pi\beta/4)}\bigg{]}.$
(71)

Here we used the normalization condition $\sum_{\tau=1}^{T}f_{\tau}=1$. The
expression for $\epsilon<0$ follows from the residue formula
$F(\epsilon)+F(-\epsilon)=2\pi$.

The ratio

$\displaystyle\bigg{|}\frac{F_{h,T}(+\epsilon)}{F_{h,T}(-\epsilon)}\bigg{|}=\bigg{|}\frac{\tanh(\epsilon\beta/2)-\tanh((\epsilon-\pi)\beta/2)}{\tanh(\epsilon\beta/2)-\tanh((\epsilon+\pi)\beta/2)}\bigg{|}$
(72) $\displaystyle\approx e^{\epsilon\beta},$

gives the detailed balance condition when $\pi\beta\gg{\rm
max}(1,\epsilon\beta)$. This condition ensures heating transitions are
exponentially less frequent than the cooling transitions of the same energy.

In Fig. 7 we compare $|F_{h,T}(\epsilon)|$ obtained from the integral
approximation, (71), and the exact sum (69). We take $\beta=10$ and show the
summation result for $T=2\beta$. Even for this relatively small value of
$T/\beta$ (or, equivalently, relatively large $\alpha=0.5$), Eq. (71) clearly
captures the main behaviour of the filter function. The main discrepancy is
due to small oscillations (‘ringing’) on top of the integral result. These
finite-time artifacts can be viewed as errors arising from the truncation of
the Fourier series limits in Eq. (69) at finite $T$. For large $T/\beta$, this
ringing error is bounded approximately as

$\displaystyle|F_{h,T}(\epsilon)-F_{h,\infty}(\epsilon)|\lesssim\frac{\exp(-\pi
T/\beta)}{A}.$ (73)

In realistic systems the errors due to noise will typically be much larger
than the above truncation errors, and we can obtain efficient cooling with
values of $T$, $\beta$, considered in the main text.

Figure 7: Comparison of integral approximation for $F_{h,T}(\epsilon)$ (Eq.
71, black) and exact summation (Eq. 69, orange), for parameters $\beta=10$,
$T=2\beta$.

##  Appendix B Detailed balance and thermal state preparation with MCP

Figure 8: Steady state quasiparticle populations for cooling with thermal MCP,
for effective temperature $T_{\rm eff}=1/2\beta$. Results shown are for PM
phase with one auxiliary per site, with system size $N_{S}=30$ and reset time
$T=5\beta$. The dashed black lines are the expected thermal (Fermi)
distributions.

Above, we focused on preparation of ground states, which corresponds to the
limit $\beta\Delta\to\infty$, where $\Delta$ is the quasiparticle gap.
However, our MCP protocol also provides a way to prepare _thermal states_
(finite temperature Gibbs states) of the same effective Hamiltonian, a subject
which has been the focus of several recent studies [27, 74, 76, 34]. The
thermal Gibbs state at an effective temperature $T_{\rm eff}$ is defined by

$\frac{\rho_{\vec{\alpha}\vec{\alpha}}}{\rho_{\vec{\beta}\vec{\beta}}}=\exp\bigg{(}-\frac{\varepsilon(\vec{\alpha})-\varepsilon(\vec{\beta})}{T_{\rm
eff}}\bigg{)},$ (74)

where $\rho_{\vec{\alpha}\vec{\alpha}}$ are the many-body populations in Eq.
(19), $\varepsilon$ is the total quasienergy of the many-body eigenstate.

This can be achieved by cooling the system using the MCP with the parameter
$\beta=1/2T_{\rm eff}$, and the additional constraint that the coupling
operators $V_{j}$ appearing in Eq. (8) are _Hermitian_. Then, the general rate
equation, Eq. (23), replacing $\hat{\sigma}^{\pm}_{j}$ with the Hermitian
operators $V_{j}$, can be written as

$\delta
n_{k}=\frac{1}{2}\sum_{\vec{\alpha}\vec{\beta}}(\beta_{k}-\alpha_{k})w_{\vec{\alpha}\vec{\beta}}\bigg{(}\rho_{\vec{\alpha}\vec{\alpha}}-\rho_{\vec{\beta}\vec{\beta}}\frac{w_{\vec{\beta},\vec{\alpha}}}{w_{\vec{\alpha}\vec{\beta}}}\bigg{)}$
(75)

where

$w_{\vec{\alpha}\vec{\beta}}=|F_{h,T}(\Delta(\vec{\alpha},\vec{\beta}))|^{2}|\bra{\vec{\beta}}V_{j}\ket{\vec{\alpha}}|^{2},$
(76)

and
$\Delta(\vec{\alpha},\vec{\beta})=\varepsilon(\vec{\alpha})-\varepsilon(\vec{\beta})$
are the Bohr frequencies. If the _detailed balance_ condition holds, that is,
when the terms inside the bracket above vanish identically,

$\frac{\rho_{\vec{\alpha}\vec{\alpha}}}{\rho_{\vec{\beta}\vec{\beta}}}=\frac{|F_{h,T}(\Delta(\vec{\beta},\vec{\alpha}))|^{2}|\bra{\vec{\alpha}}V_{j}\ket{\vec{\beta}}|^{2}}{|F_{h,T}(\Delta(\vec{\alpha},\vec{\beta}))|^{2}|\bra{\vec{\beta}}V_{j}\ket{\vec{\alpha}}|^{2}},$
(77)

then this provides a solution to the steady state $\delta n_{k}=0$. With the
Hermitian constraint and choice $\beta=1/2T_{\rm eff}$, using the asymptotic
form of the filter function Eq. (72) gives the required Gibbs state, Eq. (74).

As a demonstration, we calculate the steady state quasiparticle populations
from the rate equation in the PM phase of the Floquet TFIM, with one auxiliary
per site and keeping single-particle processes as in Section III. We choose
hermitian operators $V_{j}=\hat{Z}_{j}$. The results are shown in Fig. 8, for
several values of the inverse temperature $\beta$. We take system size
$N_{S}=30$ and $T=5\beta$ to reduce ringing artifacts. The results are
compared to the expected Fermi distributions which hold for fermionic
quasiparticles when Eq. (74) is satisfied; the agreement is generally
excellent, suggesting that thermal Gibbs states can be efficiently prepared by
the quasiparticle cooling algorithm (where, as in the ground state case,
errors due to noise etc. will result in slightly non-equilibrium steady state
populations).

##  Appendix C Quasiparticles in the TFIM

Figure 9: Profile of edge Bogoliubov coefficients $|u_{1k}|^{2}$, (solid)
$|v_{1k}|^{2}$ (dashed), plotted vs quasienergy $\epsilon_{k}$. We plot the
coefficients rescaled by chain length $N_{S}$: squared coefficients are
$O(1/N_{S})$ in middle of quasiparticle band, but vanish as $O(1/N_{S}^{3})$
at band edges.

The Floquet TFIM with open boundary conditions, considered as our cooling
example throughout the main text, is defined by the Floquet unitary

$\displaystyle\hat{U}_{S}=\exp\bigg{(}{-\frac{i\pi
J}{2}\sum_{i=1}^{N_{S}-1}\hat{X}_{i}\hat{X}_{i+1}}\bigg{)}\exp\bigg{(}\frac{i\pi
g}{2}\sum_{i=1}^{N_{S}}\hat{Z}_{i}\bigg{)}.$ (78)

This model can be mapped onto a quadratic fermionic chain by the Jordan-Wigner
mapping, Eq. (32). For this purpose it is convenient to introduce a basis of
Hermitian Majorana fermion operators, defined as

$\hat{a}_{2i-1}=\hat{c}_{i}+\hat{c}^{\dagger}_{i},\hskip
14.22636pt\hat{a}_{2i}=i(\hat{c}^{\dagger}_{i}-\hat{c}_{i}).$ (79)

Then the unitary $\hat{U}_{S}$ acts on Majorana operators as a linear
transformation:

$\hat{U}_{S}^{\dagger}\hat{a}_{i}\hat{U}_{S}=\sum_{j=1}^{2N_{S}}(K_{S})_{ij}\hat{a}_{j}.$
(80)

From Eq. (78) we find the $K_{S}$ matrix given by

$\displaystyle K_{S}=\begin{pmatrix}1&0&0&\ldots&0\\\ 0&\cos\pi J&-\sin\pi
J&\ldots&0\\\ 0&\sin\pi J&\cos\pi J&\ldots&0\\\
\vdots&\vdots&\vdots&\ddots&0\\\ 0&0&0&0&1\end{pmatrix}\begin{pmatrix}\cos\pi
g&\sin\pi g&0&0&\ldots\\\ -\sin\pi g&\cos\pi g&0&0&\ldots\\\ 0&0&\cos\pi
g&\sin\pi g&\ldots\\\ 0&0&-\sin\pi g&\cos\pi g&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{pmatrix}.$ (81)

The eigenmodes of $\hat{U}_{S}$ are Dirac fermions satisfying

$\hat{U}_{S}^{\dagger}\hat{\eta}_{k}\hat{U}_{S}=e^{-i\epsilon_{k}}\hat{\eta}_{k},$
(82)

where $k$ is a quantum number labeling the eigenmodes (the quasimomentum).
Defining the eigenvectors of $K_{S}$ (for $\epsilon_{k}\geq 0$) according to

$K_{S}\psi_{k}=e^{-i\epsilon_{k}}\psi_{k},$ (83)

we have the expansion

$\hat{\eta}_{k}=\frac{1}{\sqrt{2}}\sum_{j=1}^{N_{S}}\psi^{*}_{k,2j-1}\hat{a}_{2j-1}+\psi^{*}_{k,2j}\hat{a}_{2j}.$
(84)

Due to particle-hole symmetry, the eigenvectors of $K_{S}$ with negative
quasienergy are related to those of positive quasienergy by the conjugation
$\varphi_{k}=\psi_{k}^{*}$, $K_{S}\varphi_{k}=e^{i\epsilon_{k}}\varphi_{k}$.
The eigenvectors $\psi_{k}$ and quasienergies $\epsilon_{k}$ can be found by
exact diagonalisation of the $K_{S}$ matrix; in Appendix A of Ref. [38], an
analytical derivation was also given. Note that the boundary condition defines
the quantization condition for the $N_{S}$ quasimomenta $k_{m}$, in terms of
which the quasienergies are given by the exact formula

$\cos\epsilon_{k}=\cos(\pi J)\cos(\pi g)-\sin(\pi J)\sin(\pi g)\cos k.$ (85)

Along with the definition of the Bogoliubov coefficients,

$\hat{c}_{j}=\sum_{k}u_{jk}\hat{\eta}_{k}+v_{jk}\hat{\eta}_{k}^{\dagger},$
(86)

this gives the relations

$\displaystyle u_{jk}=\frac{\psi_{k,2j-1}+i\psi_{k,2j}}{\sqrt{2}},$ (87)
$\displaystyle v_{jk}=\frac{\psi^{*}_{k,2j-1}+i\psi^{*}_{k,2j}}{\sqrt{2}},$
(88)

specifying the standing-wave form of the eigenmodes.

Since the Bogoliubov coefficients at the boundary play a particularly
important role in the cooling rates entering the kinetic equation for edge
cooling, we plot the values of $|u_{1k}|^{2}$, $|v_{1k}^{2}|$ in Fig. 9, vs.
the quasienergy $\epsilon_{k}$. We observe that quasiparticle in the middle of
the band typically have boundary overlap $|u_{1k}|^{2}=O(N_{S}^{-1})$, while
the overlap vanishes for quasiparticles with quasienergy near the two band
edges. The overlap for these quasiparticles is expected to go as
$O(N_{S}^{-3})$ [59].

##  Appendix D Matrix elements for bulk cooling

In Sections III, IV, we considered cooling in the example of the transverse-
field Ising model, with two types of matrix elements appearing in the
derivation of the rate equations: the $\hat{Z}_{j}$ matrix elements, which we
denote as

$Z_{j;\vec{\alpha};\vec{\beta}}=\bra{\vec{\alpha}}\hat{Z}_{j}\ket{\vec{\beta}},$
(89)

and the $\hat{\sigma}^{\pm}_{j}$ matrix elements,

$\sigma^{\pm}_{j;\vec{\alpha};\vec{\beta}}=\bra{\vec{\alpha}}\hat{\sigma}^{\pm}_{j}\ket{\vec{\beta}}.$
(90)

Here $\ket{\vec{\alpha}}$ are the many-body eigenstates of the TFIM, which are
completely defined by their quasiparticle occupation numbers.

Although TFIM defines an integrable model, calculation of spin matrix elements
can be challenging [77], due to the fact that the Jordan-Wigner transformation
between spin and fermion operators is non-local. We will see that this is the
case for the $\hat{\sigma}^{\pm}_{j}$ matrix elements below. First, we focus
on the simpler case of $\hat{Z}_{j}$:

###  D.1 $Z_{j;\vec{\alpha};\vec{\beta}}$ matrix elements

In the fermion basis $\hat{Z}$ is a quadratic operator, so only matrix
elements between eigenstates differing in the occupation of at most two
quasiparticle modes survive:

$\displaystyle
Z_{j;\vec{\alpha};\vec{\beta}}=\bra{\vec{\alpha}}(1-2\hat{c}^{\dagger}_{j}\hat{c}_{j})\ket{\vec{\beta}}$
$\displaystyle=\delta_{\vec{\alpha}\vec{\beta}}-2\sum_{kq}\Big{(}u^{*}_{jk}u_{jq}[\hat{\eta}_{k}^{\dagger}\hat{\eta}_{q}]_{\vec{\alpha}\vec{\beta}}+u^{*}_{jk}v_{jq}[\hat{\eta}_{k}^{\dagger}\hat{\eta}^{\dagger}_{q}]_{\vec{\alpha}\vec{\beta}}+v^{*}_{jk}u_{jq}[\hat{\eta}_{k}\hat{\eta}_{q}]_{\vec{\alpha}\vec{\beta}}+v^{*}_{jk}v_{jq}[\hat{\eta}_{k}\hat{\eta}^{\dagger}_{q}]_{\vec{\alpha}\vec{\beta}}\Big{)}$
(91)

where e.g.
$[\hat{\eta}^{\dagger}_{k}\hat{\eta}_{q}]_{\vec{\alpha}\vec{\beta}}\equiv\bra{\vec{\alpha}}\hat{\eta}^{\dagger}_{k}\hat{\eta}_{q}\ket{\vec{\beta}}=0,1$
is only non-zero if $\vec{\alpha}$ and $\vec{\beta}$ differ only by the
occupations of the $k$ and $q$ quasiparticle levels (including the case
$k=q$). Since the Bogoliubov coefficients $u_{jk}$, $v_{jk}$ can be obtained
by diagonalisation of the single-particle Hamiltonian, the matrix elements are
easily obtained.

In deriving the rate equation, we focused on the case of transitions to and
from the ground state, and scattering between single quasiparticle levels.
Thus we have

$|Z_{j;0;k,q}|^{2}=|\bra{\Omega}\hat{Z}_{j}\ket{k,q}|^{2}=4|v_{jk}^{*}u_{jq}|^{2}$
(92)

as the probability for two quasiparticles with quasimomenta $k,\ q$, to
annihilate, and

$|Z_{j;k;q}|^{2}=|\bra{k}\hat{Z}_{j}\ket{q}|^{2}=4|u_{jk}^{*}u_{jq}|^{2}$ (93)

as the scattering probability $q\to k$. Since $\hat{Z}$ is hermitian, the
pair-creation probability is equal to the pair-annihilation, i.e.
$|Z_{j;0;k,q}|^{2}=|Z_{j;k,q;0}|^{2}$. Note that, in the AFM phase there is
also an amplitude for single-quasiparticle processes, due to the selection
rules discusses in the main text i.e. one of $k$, $q$ can be taken as the
Majorana zero mode.

###  D.2 $\sigma^{\pm}_{j;\alpha;\beta}$ matrix elements

We now turn to quasiparticle processes mediated by the
$\hat{\sigma}^{\pm}_{j}$ matrix elements. We focus on $\hat{\sigma}^{-}_{j}$
for illustration. In the fermion basis, $\hat{\sigma}^{-}_{j}$ comes attached
with the Jordan-Wigner string operator,

$\sigma^{-}_{j;\vec{\alpha};\vec{\beta}}=\bra{\vec{\alpha}}\Big{(}\prod_{i=1}^{j-1}e^{i\pi\hat{c}^{\dagger}_{i}\hat{c}_{i}}\Big{)}\hat{c}^{\dagger}_{j}\ket{\vec{\beta}}.$
(94)

Due to this non-locality in the fermionic eigenbasis, multi-particle processes
are allowed with $\hat{\sigma}_{j}^{-}$ changing the occupation numbers of up
to $2j-1$ quasiparticle levels. Amazingly, the exact finite-volume matrix
elements have been worked out for the case of the $\hat{X}_{j}$ operators in
the TFIM with periodic boundary conditions [77]. Here, we focus on singe- and
two-particle transitions near the ground state, which capture the dominant
cooling processes at late times, as described in the main text. Our theory can
be straightforwardly generalised to include higher-order processes.

Focusing on few-particle processes alows us to use Wick’s theorem to evaluate
the necessary matrix elements for systems up to several hundred sites: For
free (more generally, Gaussian) theories, Wick’s theorem states that

$\displaystyle\langle\hat{a}_{i_{1}}\hat{a}_{i_{2}}\ldots\hat{a}_{i_{2N}}\rangle_{\Omega}=\sum_{P}\zeta_{P}\prod_{k=1}^{N}\langle\hat{a}_{P(i_{2k-1})}\hat{a}_{P(i_{2k})}\rangle$
(95) $\displaystyle=\text{Pf}\
\begin{pmatrix}0&\langle\hat{a}_{i_{1}}\hat{a}_{i_{2}}\rangle&\langle\hat{a}_{i_{1}}\hat{a}_{i_{3}}\rangle&\ldots&\langle\hat{a}_{i_{1}}\hat{a}_{i_{2N}}\rangle\\\
{}&0&\langle\hat{a}_{i_{2}}\hat{a}_{i_{3}}\rangle&\ldots&\langle\hat{a}_{i_{2}}\hat{a}_{i_{2N}}\rangle\\\
{}&{}&0&\ldots&\langle\hat{a}_{i_{3}}\hat{a}_{i_{2N}}\rangle\\\
{}&{}&{}&\ddots&\vdots\\\ {}&{}&{}&{}&0\end{pmatrix}.$ (96)

The sum runs over all possible permutations of the Majorana indices, and
$\zeta_{P}=\pm 1$ is the sign of the permutation, carrying the fermion
statistics. All expectation values are taken with respect to the ground state
of the free theory, and for notational simplicity we introduced the Majorana
fermion operators (Eq. (79)). In the second line we introduced the Pfaffian
representation of Wick’s theorem in terms of a skew-symmetric matrix $\Sigma$.
In our case, since only the modulus squared of the correlation functions
enters the rate equation, we may use

$|\langle a_{i_{1}}a_{i_{2}}\ldots
a_{i_{2N}}\rangle_{\Omega}|^{2}=|\text{Pf}(\Sigma)|^{2}=|\text{det}(\Sigma)|.$
(97)

which can be computed in time $O(N^{3})$. For further reference, see e.g. [78,
61].

Single quasiparticle processes: We first consider single-quasiparticle
transitions above the ground state,

$\sigma^{-}_{j;0;k}=\bra{\Omega}\hat{\sigma}^{-}_{j}\ket{k}$ (98)

being the amplitude to directly excite a quasiparticle of quasimomentum $k$.
We rewrite the single-particle element in fermion basis

$|\sigma^{-}_{j;0;k}|^{2}=|\langle\prod_{i=1}^{j-1}\hat{a}_{2i-1}\hat{a}_{2i}\
\hat{c}^{\dagger}_{j}\hat{\eta}^{\dagger}_{k}\rangle_{\Omega}|^{2}.$ (99)

The corresponding Wick’s matrix $\Sigma$ had dimension $2j\times 2j$. In Fig.
10a we show how the single-particle matrix elements in Eq. 99 (summed over
$k$), depend on the TFIM phase parameter $J/g$ and on the position of the
operator $j$. The system size is $N_{S}=100$. The matrix elements display a
clear change in behaviour between the two phases: in the PM phase, single-
particle transitions occur almost uniformly throughout the system bulk; in the
AFM phase, transitions are localized within a few correlation lengths of the
edge, with a correlation length $\xi_{J}\propto J/g$. This is as expected,
based on the physical picture of domain-wall-like quasiparticles in the AFM
(removing a single domain wall requires flipping a macroscopic number of
spins).

Figure 10: (a) Single-quasiparticle matrix elements $|\sigma^{-}_{j0k}|^{2}$
summed over quasimomenta $k$ and plotted against operator position $j$, for
different values of $J/g$. (b) Two-particle matrix elements
$|\sigma^{-}_{j0(k,q)}|^{2}$ summed over $k$, $q$, and plotted against
operator position $j$. Note that the two-particle matrix element is non-zero
only in the AFM phase.

Two-quasiparticle processes: Based on the physical picture, we expect that in
the AFM phase, a pair of domain walls may collide in the bulk and annihilate,
as this requires only local reconfigurations of spins. However, there is a
subtlety, from the fact that $\hat{\sigma}^{-}_{j}$ has non-zero matrix
elements only between eigenstates of different fermion parity. This rules out
two-particle transitions in the PM phase, but as discussed in the main text,
the issue is sidestepped in the AFM by interacting with the Majorana zero
mode, whose action maps between the two nearly-degenerate ground states. We
therefore write

$|\sigma^{-}_{j;0;k,q}|^{2}=|\bra{\Omega_{+}}\prod_{i=1}^{j-1}\hat{a}_{2i-1}\hat{a}_{2i}\
\hat{c}^{\dagger}_{j}\hat{\eta}^{\dagger}_{k}\hat{\eta}^{\dagger}_{q}\ket{\Omega_{-}}|^{2},$
(100)

with $\ket{\Omega_{-}}=\hat{\eta}^{\dagger}_{0}\ket{\Omega_{-}}$, and
$\hat{\eta}^{\dagger}_{0}$ the Majorana mode. These matrix elements may be
evaluated with Wick’s theorem on the $(2j+2)\times(2j+2)$ dimensional
$\Sigma$-matrix. In Fig. 10b, we show how the summed matrix elements
$\sum_{k,q}|\sigma^{-}_{j;0;k,q}|^{2}$ depend on the position $j$, for
different values of $J/g$ in the AFM phase. The pair matrix element is almost
homogeneous in the bulk, vanishing in the edge region.

For completeness we also give the two-particle scattering elements, which are
again non-zero only for the AFM phase:

$|\sigma^{-}_{j;k;q}|^{2}=|\bra{\Omega_{+}}\hat{\eta}_{k}\prod_{i=1}^{j-1}\hat{a}_{2i-1}\hat{a}_{2i}\
\hat{c}^{\dagger}_{j}\hat{\eta}^{\dagger}_{q}\ket{\Omega_{-}}|^{2}.$ (101)

The scattering does not directly lead to cooling/heating of the system (the
total number of quasiparticles is conserved by the scattering processes),
however it can lead to a redistribution of quasiparticle density, which itself
affects the cooling rate.

##  Appendix E Details of numerical simulations

###  E.1 Free fermion simulations

For the edge cooling protocol analysed in Section III for the TFIM, the non-
unitary dynamics described by Eq. (1) defines a Gaussian channel — a channel
mapping Gaussian states to Gaussian states. A Gaussian state obeys Wick’s
theorem (see Eq. 96), meaning that the expectation values of all observables
can be factorised purely in terms of $2$-point fermionic correlation
functions. Under the action of the cooling channel, Eq. (1), the matrix of
2-point functions can be efficiently updated, as we show below, and hence we
have access to the arbitrary correlators of the time-dependent many-body
state.

We introduce the antisymmetrised covariance matrix of Majorana correlators as

$\Gamma_{ij}=\frac{i}{4}\langle[\hat{a}_{i},\hat{a}_{j}]\rangle=\begin{cases}i\langle\hat{a}_{i}\hat{a}_{j}\rangle/2,&i\neq
j\\\ 0,&i=j\end{cases},$ (102)

where the Majorana fermions were defined in Eq. (79). The matrix includes both
system and bath degrees of freedom, and has dimension
$(2N_{S}+2)\times(2N_{S}+2)$ (we limit our discussion to the case of a single
edge auxiliary for simplicity). We block the covariance matrix according to

$\Gamma=\begin{pmatrix}\Gamma_{BB}&\Gamma_{SB}\\\
\Gamma_{BS}&\Gamma_{SS}\end{pmatrix}.$ (103)

where $\Gamma_{SS}$ is the $2N_{S}\times 2N_{S}$ matrix containing system-
system correlations, etc. For the initial condition, representing the system
in the maximally mixed state and the auxiliary in the reset state, Eq. (6), we
have

$\Gamma_{SS}=0$ (104) $\Gamma_{BB}=\frac{1}{2}\begin{pmatrix}0&-1\\\
1&0\end{pmatrix},\hskip 14.22636pt\Gamma_{SB}=\Gamma_{BS}=0.$ (105)

Due to the integrability of the TFIM, unitary evolution of the Majoranas in
the Heisenberg picture acts as a linear map:

$(\hat{U}^{\dagger}_{S}\hat{U}^{\dagger}_{B}\hat{U}^{\dagger}_{\theta})\
\hat{a}_{i}\
(\hat{U}_{\theta}\hat{U}_{B}\hat{U}_{S})=\sum_{j}(K_{\theta}K_{B}K_{S})_{ij}\hat{a}_{j},$
(106)

where $K_{S}$ is defined in (81), and $K_{B}=\exp(\pi h_{B})$,
$K_{\theta}=\exp(\pi h_{\theta})$, with

$h_{B}=\begin{pmatrix}0&h&0&\ldots\\\ -h&0&0&\ldots\\\ 0&0&0&\ldots\\\
\vdots&\vdots&\vdots&\ddots\end{pmatrix},\hskip
8.5359pth_{\theta}=\begin{pmatrix}0&0&0&\theta&\ldots\\\
0&0&-\theta&0&\ldots\\\ 0&\theta&0&0&\ldots\\\ -\theta&0&0&0&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{pmatrix},$ (107)

all other elements being zero. Then, defining
$\hat{U}=\hat{U}_{\theta}\hat{U}_{B}\hat{U}_{S}$ and $K=K_{\theta}K_{B}K_{S}$,
the covariance matrix evolves under the unitary part of the cooling channel
as,

$\langle\hat{U}^{\dagger}\hat{\Gamma}_{ij}\hat{U}\rangle=K_{il}K_{jm}\langle[\hat{a}_{l},\hat{a}_{m}]\rangle=(K\Gamma
K^{T})_{ij}.$ (108)

Figure 11: Tensor network implementation for one cooling cycle. System
(auxiliary) qubits are denoted by black (red) tensors respectively. Different
channels correspond to different blocks: bath evolution (light grey); system
evolution (dark blue); system-bath coupling (dark grey); decoherence (yellow);
auxiliary reset (white). The network is composed of rank-3 tensors $A$, has a
bond dimension $d$ and a 4-dimensional physical bond corresponding to the four
(vectorised) states of the qubit’s density matrix. The system is encoded in a
one-dimensional geometry at the price of introducing next-nearest neighbor
channels.

The reset of the auxiliaries can be done by reinstating condition (105). If we
denote the corresponding map by $\mathcal{R}$, then the cooling channel (1)
acts on the covariance matrix according to (in a slight abuse of notation):

$[\Phi(\Gamma)]_{ij}=[\mathcal{R}(K\Gamma K^{T})]_{ij}.$ (109)

Thus the matrix $\Gamma$ can be efficiently updated to obtain the dynamics, or
(by vectorising the matrix $\Gamma$) the linear map $\Phi$ diagonalised to
yield the steady state. Note that qubit dephasing can also be efficiently
implemented in this scheme [79], but bulk auxiliaries or qubit decay break the
necessary integrability.

###  E.2 Tensor network simulations

Here, we summarize the tensor network algorithm that was employed to simulate
the dissipative dynamics in Sections III-V. We give an overview of the
relevant details, assuming the reader is familiar with the basics of matrix
product algorithms; introductions can be found in e.g. [80, 81].

We assume a finite density of auxiliaries, one per system site. The total
density matrix of the system plus auxiliaries is $\rho_{SB}$. We work with the
vector unfolding of the density matrix product operator, which is written as a
matrix product state (MPS),

$\rho_{SB}=\sum_{\textbf{m},\textbf{s}}A^{1}_{m_{1},s_{1}}A^{2}_{m_{2},s_{1},s_{2}}\ldots
A^{N}_{m_{N},s_{N-1}}|m_{1}\ldots m_{N}\rangle,$ (110)

where $A^{i}$ is a rank-3 tensor graphically defined in Fig. 11. Virtual
indices $\textbf{s}=\\{s_{1},\ldots,s_{N-1}\\}$ are traced while physical
indices $\textbf{m}=\\{m_{1},\ldots,m_{N}\\}$ denote the different states of
the system+bath, $N=2N_{S}$. The system geometry in the presence of
auxiliaries is quasi-1d, see Fig. 1d, and we adopt a 1-dimensional ordering
that alternates between system and auxiliary qubits i.e.
$\mathcal{Q}_{1}\mathcal{A}_{1}\ldots\mathcal{Q}_{N_{S}}\mathcal{A}_{N_{S}}$,
as illustrated in Fig. 11. The dimension of the virtual indices of the MPS
(bond dimension) is denoted $d$, and controls the amount of non-local
correlations that can be represented by the state.

Figure 12: Convergence of ground state cooling simulations as a function of
the MPS’s bond dimension, for the integrable TFIM with the parameters used in
Figure 3. We show data for both PM and AFM phases at two different bath
couplings $\theta=0.1$ $(0.05)$, and for bond dimensions $d=100,150,200,250$.
Figure 13: Convergence of ground state cooling simulations as a function of
the MPS’s bond dimension, for the non-integrable Ising model with the
parameters used in Figure 6. Full/dashed/dotted bond dimensions
$d=200,300,400$ respectively. For increasing system sizes, there is a small
drift with $d$ at $\nu>10$, but the monotonic increase of fidelity is observed
for all bond dimensions.

The time evolution is performed via the standard TEBD algorithm [64, 65],
which updates the MPS in real time, followed by the reset channel acting on
the bath, which can be efficiently implemented as discussed below. In Fig. 11
we illustrate the evolution for the example of the Floquet TFIM (5) with bath
evolution and interaction, Eqs. (7-10). We first simultaneously perform the
single-site unitary channels for the system and bath, then the two-body system
channels, which correspond to next-nearest-neighbor gates, followed by the
nearest-neighbor system-bath channel. Note that the next-nearest-neighbour
gate arises from the quasi-1d geometry. This is the most expensive step of the
algorithm as it requires two singular value decompositions and scales as
$O(\text{dim}(m)^{4}d^{3})$.

For the simulations with decoherence, following the unitary channel we apply
the decoherence channel Eq. (49) with jump operators in Eq. (54). Finally the
reset of auxiliaries is performed by simply modifying the matrix elements of
each auxiliary qubit tensor as,

$\displaystyle A^{i}_{00s_{i-1}s_{i}}\rightarrow
A^{i}_{00s_{i-1}s_{i}}+A^{i}_{11s_{i-1}s_{i}},$ (111) $\displaystyle
A^{i}_{m\neq 00s_{i-1}s_{i}}\rightarrow 0.$

In the non-integrable TFIM in the presence of a longitudinal field, the local
unitary operations are defined according to Eq. (66), and the rest of the
algorithm is unaffected.

Next, we study the effects of truncation due to finite bond dimension during
the cooling protocol in Figures 12 and 13. In Fig. 12 we consider the TFIM
cooling for the parameters shown in Fig. 3 of the main text. We only show data
for two values of the coupling, $\theta=0.1,0.05$ for the MCP, since the SCP
protocol is well converged for all bond dimensions. We observe that the
convergence is slower for the larger coupling, suggesting higher bond
dimensions are needed to capture non-perturbative effects. Phenomenologically,
such behavior is expected since higher-order in $\theta$ corrections are
expected to lead processes that are less spatially local, thus requiring a
larger bond dimension to be accurately captured.

The comparison of data obtained with different bond dimension $d=200,300,400$
for the non-integrable model is shown in Fig. 13. We note that the bond
dimension necessary for convergence depends on the correlation length of the
ground state. For the parameters of the TFIM with longitudinal field used in
the main text ($h_{x}=h_{z}=0.15,J=0.1$), the effects of truncation are almost
negligible for bond dimension values considered.

##  Appendix F Comparison of noisy kinetic equation and MPS simulations

In Section IV we derived the kinetic equation in the presence of noise, Eq.
(51), which includes transitions induced by qubit dephasing and decay, Eqs.
(51). Here we compare the accuracy of the rate equation predictions against
MPS simulations. We focus on the MCP with the same parameters as in Fig. (4),
namely finite density of auxiliaries, $N_{A}=N_{S}=20$, and $T=12=2\beta$. In
Fig. 14 we plot the quasiparticle densities in the steady state, obtained from
the kinetic equation prediction (black dashed line) and MPS simulations
(scatter points). We use a maximum MPS bond dimension of $d=100$ and run the
simulations starting from the exact ground state, as described in Section III.
We fix the noise strength $\gamma_{d}=\gamma_{\phi}=0.1\theta^{2}$, and show
results for three different coupling strengths, $\theta=0.2,0.1,0.05$, finding
good qualitative agreement for the smaller coupling values $\theta=0.1,0.05$.

Figure 14: Quasiparticle occupation numbers $n^{\infty}_{k}$ in steady state,
as a function of quasienergy $\epsilon_{k}$, for noise strength
$\gamma_{d}=\gamma_{\phi}=0.1\theta^{2}$. (a) PM and (b) AFM phases. Rate
equation prediction (dashed line) vs. MPS simulations (points), for three
values of the coupling $\theta$. Protocol parameters are as in Fig. (4). The
MPS simulations use a maximum bond dimension of $d=100$.

## References

  * Lanyon  _et al._ [2010] B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, _et al._ , Towards quantum chemistry on a quantum computer, Nature chemistry 2, 106 (2010).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, nature 549, 242 (2017).
  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Quantum algorithms for quantum chemistry and quantum materials science, Chemical Reviews 120, 12685 (2020).
  * Quantum  _et al._ [2020] G. A. Quantum, Collaborators*†, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, _et al._ , Hartree-fock on a superconducting qubit quantum computer, Science 369, 1084 (2020).
  * Finnila  _et al._ [1994] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll, Quantum annealing: A new method for minimizing multidimensional functions, Chemical physics letters 219, 343 (1994).
  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem, Science 292, 472 (2001).
  * Kitaev [1995] A. Y. Kitaev, Quantum measurements and the abelian stabilizer problem (1995), arXiv:quant-ph/9511026 [quant-ph] .
  * Abrams and Lloyd [1999] D. S. Abrams and S. Lloyd, Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors, Phys. Rev. Lett. 83, 5162 (1999).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum computation by adiabatic evolution, arXiv preprint quant-ph/0001106 (2000).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Preparing ground states of quantum many-body systems on a quantum computer, Phys. Rev. Lett. 102, 130503 (2009).
  * Boixo  _et al._ [2009] S. Boixo, E. Knill, and R. D. Somma, Eigenpath traversal by phase randomization, Quantum Inf. Comput. 9, 833 (2009).
  * Ge  _et al._ [2019] Y. Ge, J. Tura, and J. I. Cirac, Faster ground state preparation and high-precision ground energy estimation with fewer qubits, Journal of Mathematical Physics 60, 022202 (2019).
  * Preskill [2018] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Reviews of Modern Physics 90, 015002 (2018).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * Simon  _et al._ [2011] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Quantum simulation of antiferromagnetic spin chains in an optical lattice, Nature 472, 307 (2011).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021), https://www.science.org/doi/pdf/10.1126/science.abi8794 .
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms, Nature 595, 233 (2021).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The variational quantum eigensolver: A review of methods and best practices, Physics Reports 986, 1 (2022).
  * Fedorov  _et al._ [2022] D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, Vqe method: a short survey and recent developments, Materials Theory 6, 1 (2022).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Louvet  _et al._ [2023] T. Louvet, T. Ayral, and X. Waintal, Go-no go criteria for performing quantum chemistry calculations on quantum computers (2023), arXiv:2306.02620 [quant-ph] .
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature communications 9, 4812 (2018).
  * Wang  _et al._ [2021] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise-induced barren plateaus in variational quantum algorithms, Nature Communications 12, 6961 (2021).
  * Sack  _et al._ [2022] S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, Avoiding barren plateaus using classical shadows, PRX Quantum 3, 020365 (2022).
  * Terhal and DiVincenzo [2000] B. M. Terhal and D. P. DiVincenzo, Problem of equilibration and the computation of correlation functions on a quantum computer, Phys. Rev. A 61, 022301 (2000).
  * Davies [1974] E. B. Davies, Markovian master equations, Communications in mathematical Physics 39, 91 (1974).
  * Kraus  _et al._ [2008] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Preparation of entangled states by quantum markov processes, Phys. Rev. A 78, 042307 (2008).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Quantum states and phases in driven open quantum systems with cold atoms, Nat. Phys. 4, 878 (2008).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Quantum computation and quantum-state engineering driven by dissipation, Nat. Phys. 5, 633 (2009).
  * Roy  _et al._ [2020] S. Roy, J. Chalker, I. Gornyi, and Y. Gefen, Measurement-induced steering of quantum systems, Physical Review Research 2, 033347 (2020).
  * Kastoryano and Brandao [2016] M. J. Kastoryano and F. G. Brandao, Quantum gibbs samplers: The commuting case, Communications in Mathematical Physics 344, 915 (2016).
  * Chen  _et al._ [2023] C.-F. Chen, M. J. Kastoryano, F. G. S. L. Brandão, and A. Gilyén, Quantum thermal state preparation (2023), arXiv:2303.18224 [quant-ph] .
  * Barreiro  _et al._ [2011] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, An open-system quantum simulator with trapped ions, Nature 470, 486 (2011).
  * Lin  _et al._ [2013] Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler, A. S. Sørensen, D. Leibfried, and D. J. Wineland, Dissipative production of a maximally entangled steady state of two quantum bits, Nature 504, 415 (2013).
  * Shankar  _et al._ [2013] S. Shankar, M. Hatridge, Z. Leghtas, K. M. Sliwa, A. Narla, U. Vool, S. M. Girvin, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Autonomously stabilized entanglement between two superconducting quantum bits, Nature 504, 419 (2013).
  * Mi  _et al._ [2024] X. Mi, A. Michailidis, S. Shabani, K. Miao, P. Klimov, J. Lloyd, E. Rosenberg, R. Acharya, I. Aleiner, T. Andersen, _et al._ , Stable quantum-correlated many-body states through engineered dissipation, Science 383, 1332 (2024).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, Effective Hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Physical Review B 95, 014112 (2017).
  * Metcalf  _et al._ [2020] M. Metcalf, J. E. Moussa, W. A. de Jong, and M. Sarovar, Engineered thermalization and cooling of quantum many-body systems, Physical Review Research 2, 023214 (2020).
  * Raghunandan  _et al._ [2020] M. Raghunandan, F. Wolf, C. Ospelkaus, P. O. Schmidt, and H. Weimer, Initialization of quantum simulators by sympathetic cooling, Sci. Adv. 6, eaaw9268 (2020).
  * Polla  _et al._ [2021] S. Polla, Y. Herasymenko, and T. E. O’Brien, Quantum digital cooling, Phys. Rev. A 104, 012414 (2021).
  * Matthies  _et al._ [2023] A. Matthies, M. Rudner, A. Rosch, and E. Berg, Programmable adiabatic demagnetization for systems with trivial and topological excitations (2023), arXiv:2210.17256 [quant-ph] .
  * Puente  _et al._ [2023] D. A. Puente, F. Motzoi, T. Calarco, G. Morigi, and M. Rizzi, Quantum state preparation via engineered ancilla resetting, arXiv preprint arXiv:2305.08641 (2023).
  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * Sachdev [1999] S. Sachdev, Quantum phase transitions, Physics world 12, 33 (1999).
  * Ovchinnikov  _et al._ [2003] A. A. Ovchinnikov, D. V. Dmitriev, V. Y. Krivnov, and V. O. Cheranovskii, Antiferromagnetic ising chain in a mixed transverse and longitudinal magnetic field, Phys. Rev. B 68, 214406 (2003).
  * Note [1] For protocols considered below, $f_{\tau}$ can be positive or negative, but for this choice the sum $\sum@\slimits@_{\tau=1}^{T}|f_{\tau}|\gtrsim 1$.
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, USA, 2002).
  * Lidar [2020] D. A. Lidar, Lecture notes on the theory of open quantum systems (2020), arXiv:1902.00967 [quant-ph] .
  * Note [2] We note that the second assumption is practically necessary, as the diagonal assumption still leaves us with $O(\mathop{exp}\nolimits(N_{S}))$ coupled equations to solve.
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons, Phys. Rev. Lett. 98, 050405 (2007).
  * Lange  _et al._ [2018] F. Lange, Z. Lenarčič, and A. Rosch, Time-dependent generalized gibbs ensembles in open quantum systems, Physical Review B 97, 165138 (2018).
  * [55] We note that the normalization factor is chosen according to the condition $\sum@\slimits@_{\tau=1}^{T}f_{\tau}=1$. The modulated coupling in Eq. (29) is sign-changing, however, $\sum@\slimits@_{\tau=1}^{T}|f_{\tau}|$ is of order 1 for all cases considered below.
  * Thakurathi  _et al._ [2013] M. Thakurathi, A. A. Patel, D. Sen, and A. Dutta, Floquet generation of majorana end modes and topological invariants, Phys. Rev. B 88, 155133 (2013).
  * Bauer  _et al._ [2019] B. Bauer, T. Pereg-Barnea, T. Karzig, M.-T. Rieder, G. Refael, E. Berg, and Y. Oreg, Topologically protected braiding in a single wire using floquet majorana modes, Phys. Rev. B 100, 041102 (2019).
  * Lerose  _et al._ [2021] A. Lerose, M. Sonner, and D. A. Abanin, Scaling of temporal entanglement in proximity to integrability, Phys. Rev. B 104, 035137 (2021).
  * Žnidarič [2015] M. Žnidarič, Relaxation times of dissipative many-body quantum systems, Physical Review E 92, 042143 (2015).
  * Note [3] Heating processes, in particular due to noise, may in principle involve a large number of particles; however, we will find that few-quasiparticle processes are dominant, and considering them gives an accurate approximation of heating rates.
  * Mbeng  _et al._ [2020] G. B. Mbeng, A. Russomanno, and G. E. Santoro, The quantum ising chain for beginners (2020), arXiv:2009.09208 [quant-ph] .
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004).
  * Zwolak and Vidal [2004] M. Zwolak and G. Vidal, Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm, Phys. Rev. Lett. 93, 207205 (2004).
  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Physical review letters 91, 147902 (2003).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Physical review letters 93, 040502 (2004).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Note [4] The steady state can be efficiently found by root-finding algorithms, after the cooling and noise-induced rates have been calculated via Wick’s theorem. In practice, the system reaches its steady state quickly and we evolve the rate equation in time until the quasiparticle occupations have effectively converged.
  * Moll  _et al._ [2018] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, A. Kandala, A. Mezzacapo, P. Müller, W. Riess, G. Salis, J. Smolin, I. Tavernelli, and K. Temme, Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Modugno  _et al._ [2001] G. Modugno, G. Ferrari, G. Roati, R. J. Brecha, A. Simoni, and M. Inguscio, Bose-einstein condensation of potassium atoms by sympathetic cooling, Science 294, 1320 (2001), https://www.science.org/doi/pdf/10.1126/science.1066687 .
  * Kempe  _et al._ [2006] J. Kempe, A. Kitaev, and O. Regev, The complexity of the local hamiltonian problem, Siam journal on computing 35, 1070 (2006).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-robinson bounds and the generation of correlations and topological quantum order, Physical review letters 97, 050401 (2006).
  * Hastings [2011] M. B. Hastings, Topological order at nonzero temperature, Physical review letters 107, 210501 (2011).
  * Kishony  _et al._ [2023] G. Kishony, M. S. Rudner, A. Rosch, and E. Berg, Gauged cooling of topological excitations and emergent fermions on quantum simulators, arXiv preprint arXiv:2310.16082 (2023).
  * Chen and Brandao [2021] C.-F. Chen and F. G. Brandao, Fast thermalization from the eigenstate thermalization hypothesis, arXiv preprint arXiv:2112.07646 (2021).
  * Shtanko and Movassagh [2023] O. Shtanko and R. Movassagh, Preparing thermal states on noiseless and noisy programmable quantum processors (2023), arXiv:2112.14688 [quant-ph] .
  * Rall  _et al._ [2023] P. Rall, C. Wang, and P. Wocjan, Thermal State Preparation via Rounding Promises, Quantum 7, 1132 (2023).
  * Iorgov  _et al._ [2011] N. Iorgov, V. Shadura, and Y. Tykhyy, Spin operator matrix elements in the quantum ising chain: fermion approach, Journal of Statistical Mechanics: Theory and Experiment 2011, P02028 (2011).
  * Barouch and McCoy [1971] E. Barouch and B. M. McCoy, Statistical mechanics of the x y model. ii. spin-correlation functions, Physical Review A 3, 786 (1971).
  * Barthel and Zhang [2022] T. Barthel and Y. Zhang, Solving quasi-free and quadratic lindblad master equations for open fermionic and bosonic systems, Journal of Statistical Mechanics: Theory and Experiment 2022, 113101 (2022).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of physics 349, 117 (2014).
  * Bridgeman and Chubb [2017] J. C. Bridgeman and C. T. Chubb, Hand-waving and interpretive dance: an introductory course on tensor networks, Journal of physics A: Mathematical and theoretical 50, 223001 (2017).