Efficiently preparing chiral states via fermionic cooling on bosonic quantum hardware
Gilad Kishony Mark S. Rudner Erez Berg











    

Department of Condensed Matter Physics,
Weizmann Institute of Science,
Rehovot 76100, Israel

Department of Physics, University of Washington, Seattle, WA 98195-1560, USA

Department of Condensed Matter Physics,
Weizmann Institute of Science,
Rehovot 76100, Israel

    
§ ABSTRACT 

We propose an efficient protocol for preparing low energy states of arbitrary fermionic Hamiltonians on a noisy bosonic quantum simulator. This procedure involves performing adiabatic cooling by coupling the target system with a simulated bath. The bath is periodically monitored in order to extract entropy from the system. By fermionizing the simulated target system and the bath together, we allow individual fermionic excitations of the system to coherently hop to the bath sites. In this way, we achieve a cooling rate linearly proportional to the density of these excitations, despite the fact that they are non-local in terms of the bosonic degrees of freedom of the hardware. In particular, we show that certain topological phases, such as the chiral (non-Abelian) phase of the Kitaev honeycomb model can be prepared efficiently using our protocol. We find that our protocol performs favorably in the presence of noise, making it suitable for execution on near-term quantum devices.

   
    
Efficiently preparing chiral states via fermionic cooling on bosonic quantum hardware
    Erez Berg
    

=====================================================================================





§ INTRODUCTION



Simulating quantum many-body systems, both in and out of equilibrium, is one of most promising applications of near-term quantum computers <cit.>. In such tasks, quantum computers offer a vast advantage over their classical counterparts, by virtue of their ability to store the many-body wave functions and perform unitary gates 
that implement physical time evolution according to arbitrary Hamiltonians, at a polynomial cost in the system size. Such simulations may have practical utility in many fields, such as materials science, high-energy physics, and quantum chemistry. Preparation of a low energy state of a many-body Hamiltonian is an important subroutine for these simulations. Doing this efficiently for arbitrary Hamiltonians is a highly non-trivial task, and is one of the major challenges in the field.

Of particular interest is the case where the ground state of the Hamiltonian is topologically non-trivial. Predicting the occurrence of topological states, e.g., in quantum spin systems or interacting electrons
in flat Chern bands, is notoriously difficult. In addition, such topological states are particularly challenging to prepare on quantum computers using conventional methods since, by definition, they cannot be connected to a product state by a finite-depth spatially-local unitary circuit <cit.>. These unitary circuit algroithms include variational methods <cit.>, adiabatic processes <cit.>, and effective imaginary time evolution <cit.> . In addition, protocols for cooling the system by coupling it to a simulated thermal bath have been proposed <cit.>. All these methods are expected to 
suffer a parametric reduction in performance

when preparing topological states, requiring
the circuit depth to increase with the system size to achieve a given accuracy.



A large class
of topologically non-trivial states can be prepared efficiently using dynamic circuits with
measurements and classical feedback <cit.>. The principle is illustrated by
the example of the toric code model <cit.>: the Hamiltonian is “frustration–free,”
i.e., it consists of a sum of commuting terms that can be measured
simultaneously. After the measurement, the system can be brought into
its ground state by applying a unitary transformation designed to
remove the excitations that were detected by the measurement. This
method is not applicable for generic Hamiltonians away from the frustration-free
limit. 
Moreover, certain topological
states, such as chiral states,  


do not have a frustration-free description. 



















In this study, we tackle the problem of preparing chiral topological states, focusing on the well-known example of the chiral gapped phase of
the Kitaev spin liquid (KSL) model on the honeycomb lattice, which hosts non-Abelian excitations <cit.>. 

The KSL model can be mapped into a system of emergent Majorana fermions coupled to a static ℤ_2 gauge field. 

The preparation scheme consists of two steps. 
In the first step, the system is prepared in the flux-free sector by measuring the flux operators through all plaquettes, and applying unitary operators that correct the fluxes that are detected <cit.>.
The second step is to prepare the fermionic ground state, whose Chern number is non-zero. The key idea is to introduce auxiliary “bath” degrees of freedom that simulate an effective fermionic bath,  coupled to the system's emergent fermions. 
Both species of fermions are coupled to the same ℤ_2 gauge field, allowing single-fermion hopping between them. 
We provide an explicit encoding of the system and bath fermions into the physical bosonic qubits, generalizing the one-dimensional construction of Ref. <cit.> to two dimensions.  
The bath fermions are then used to extract energy and entropy from the system, by performing simulated cyclic adiabatic cooling <cit.> followed by measurements of the auxiliary qubits. The steps of the protocol are illustrated in Fig. <ref>.

Our protocol prepares the ground state of the chiral phase parametrically more efficiently than other methods, such as adiabatic preparation <cit.> or `naive' adiabatic cooling using a simulated bosonic bath <cit.>. These methods require a quantum circuit whose depth grows at least polynomially with the system size. In contrast, within our method, the energy density decreases exponentially with the number of cycles performed, and thus the total depth required to achieve a given accuracy in the total ground state energy grows only logarithmically with the system size. In the presence of decoherence, our protocol reaches a steady state whose energy density is proportional to the noise rate, which is parametrically lower than bosonic cooling protocols (where the energy density of the steady state scales as the square root of the noise rate <cit.>). 









The protocol presented here can be used to efficiently prepare the ground state of a generic (interacting) fermionic model in 2D. This is done by artificially introducing a ℤ_2 gauge field similar to the physical gauge field of the KSL model in order to facilitate the fermionization scheme for the system and the bath. The presence of an effective fermionic bath allows to remove single fermion excitations, while keeping the coupling between the system and auxiliary qubits spatially local.

This paper is organized as follows. In Sec. <ref> we introduce our algorithm for preparing the ground state of the KSL model. 

In Sec. <ref> we analyze the performance of the protocol both with and without noise via numerical simulations.
Finally, in Sec. <ref> we 
generalize  our protocol for arbitrary interacting fermionic systems.



§ PREPARING A CHIRAL SPIN LIQUID





 §.§ Model


In this section we explain how 
to prepare the ground state of the KSL model in the phase where the emergent fermions carry a non-trivial Chern number, 
as efficiently as in a topologically trivial phase. The same principle can be applied to arbitrary interacting fermionic Hamiltonians, as we discuss below.

The KSL model consists of spin-1/2 degrees of freedom on a honeycomb lattice with strongly anisotropic exchange interactions <cit.>. The Hamiltonian is written as:

    H_KSL=   -∑_α∈{x,y,z}𝒥_α∑_⟨𝐈,𝐉⟩∈α-bondsσ_𝐈^ασ_𝐉^α
       -κ∑_⟨⟨𝐉,𝐊,𝐋⟩⟩σ_𝐉^xσ_𝐊^yσ_𝐋^z.

Here, 𝐉 = (i,j,s) denotes the sites of the honeycomb lattice, where i ∈{1,…,N_x} and j ∈{1,…,N_y} label the unit cell and s∈{A,B} labels the two sublattices, 
and σ_𝐉^α denotes the α=x,y,z Pauli matrix at site 𝐉. 


The bonds of the lattice are partitioned into three sets – x, y, and z bonds –according to their geometric orientations (see Fig. <ref>).
We denote the corresponding Kitaev anisotropic exchange couplings by 𝒥_x, 𝒥_y, and 𝒥_z, respectively. 
The κ term is a three-spin interaction acting on three neighboring sites denoted by ⟨⟨𝐉,𝐊,𝐋⟩⟩, 
such that both 𝐉 and 𝐋 are nearest neighbors of site 𝐊. 



The κ term breaks time reversal symmetry, and opens a gap in the bulk 
that drives the system into the chiral non-Abelian phase <cit.>. 






In order to facilitate an efficient protocol to prepare the ground state of (<ref>), we introduce a modified KSL Hamiltonian, H_KSL, that includes both a system (σ) and a auxiliary (τ) spin at every site of the honeycomb lattice. The auxiliary spins are used to cool the system down to its ground state. 






The Hamiltonian is illustrated in Fig. <ref>, and is composed of a sum of two terms:

    H(t)=H_KSL+H_ c(t),

where H_KSL is the modified KSL Hamiltonian:

    H_KSL=   -∑_α∈{x,y,z}𝒥_α∑_⟨𝐈,𝐉⟩∈α-bondsσ_𝐈^ασ_𝐉^ατ_𝐈^zτ_𝐉^z
       -κ∑_⟨⟨𝐉,𝐊,𝐋⟩⟩σ_𝐉^xτ_𝐉^zσ_𝐊^yτ_𝐊^zσ_𝐋^z.

The τ^z_ J operators all commute with H_KSL. Clearly, H_KSL is identical to  H_KSL in Eq. (<ref>) in the sector τ^z_𝐉=+1. In fact, the two Hamiltonians are related by a unitary transformation for any state of the τ spins with even parity, ∏_𝐉τ^z_𝐉=1. 
To see this, 
we note that 
all of the unitary operators 
𝒰_𝐈,𝐉 = τ^x_𝐈σ^α_𝐈σ^α_𝐉τ^x_𝐉, for ⟨𝐈,𝐉⟩∈α-bonds, commute with H_KSL and anticommute individually with τ^z_𝐈 and τ^z_𝐉. The transformation from a sector with a given {τ^z_𝐉} of even parity to the sector {τ^z_𝐉=1} can be constructed by annihilating pairs of τ^z_𝐊,τ^z_𝐋=-1 using a product of the unitaries 𝒰_𝐈,𝐉 along a path connecting the points 𝐊, 𝐋. Thus, preparing the ground state of H_KSL in any known even-parity sector of {τ^z_𝐉} is equivalent to preparing the ground state of the original KSL Hamiltonian, Eq. (<ref>). 

The “control” Hamiltonian H_ c(t), which is used to cool into the ground state of H_KSL, is chosen as a time-dependent effective Zeeman field acting on τ_𝐉 with x̂ and ẑ components g_(t) and B_𝐉(t), respectively, 
 
    H_ c(t)=-∑_𝐉B_𝐉(t)τ_𝐉^z-∑_𝐉g_(t)τ_𝐉^x.

The time dependence of B_𝐉(t) and g_(t) and the protocol used for cooling are described below.




 §.§ Fermionization


It is useful to fermionize the spins σ and τ by introducing a set of six Majorana operators {c_𝐉^x,c_𝐉^y,c_𝐉^z,b_𝐉^x,b_𝐉^y,b_𝐉^z} per site, subject to the constraint

    c_𝐉^xc_𝐉^yc_𝐉^zb_𝐉^xb_𝐉^yb_𝐉^z=i.


The spin operators are related to the fermionic ones by: 


    σ_𝐉^α  =-i/2∑_β,γε^αβγb_𝐉^βb_𝐉^γ,  τ_𝐉^α  =-i/2∑_β,γε^αβγc_𝐉^βc_𝐉^γ,

where α,β,γ∈{x,y,z} and ε^αβγ is the totally antisymmetric
tensor. This transformation of the spins to fermions is analogous to the one used by Kitaev <cit.>. Note, however, that we have fermionized the σ and τ spins together, and as a result, there is a single constraint per site containing both a σ and a τ spin. 
For later use, we note that, using Eqs. (<ref>) and (<ref>), we can write

    τ_𝐉^zσ_𝐉^α=-ic_𝐉^zb_𝐉^α.


Using this mapping, we express 
Hamiltonian (<ref>) as 






    H(t)=   ∑_α∈{x,y,z}𝒥_α∑_⟨𝐈,𝐉⟩∈α-bondsu_𝐈,𝐉ic_𝐈^zc_𝐉^z
       -iκ∑_⟨⟨𝐉,𝐊,𝐋⟩⟩u_𝐉,𝐋u_𝐋,𝐊c_𝐉^zc_𝐊^z
       -∑_𝐉B_𝐉(t)ic_𝐉^yc_𝐉^x-∑_𝐉g_(t)ic_𝐉^zc_𝐉^y,

where we identify the set of conserved ℤ_2 gauge fields

    {u_𝐈,𝐉=ib_𝐈^αb_𝐉^α|⟨𝐈,𝐉⟩∈α-bonds,α∈{x,y,z}}.


The gauge-invariant flux W_i,j of the ℤ_2 gauge field is defined on each hexagonal plaquette that corresponds to the unit cell (i,j), 

as the product of u_𝐈,𝐉 on its edges. In terms of the spin degrees of freedom, this can be written as

    W_i,j = σ^x_(i,j,B)σ^z_(i+1,j,A)σ^y_(i+1,j,B)
    ×σ^x_(i+1,j+1,A)σ^z_(i,j+1,B)σ^y_(i,j+1,A).


The fermionic Hamiltonian (<ref>) with g_=0, B_𝐉=0 is identical to the Hamiltonian found by fermionizing the 
original KSL model <cit.>, with {c_𝐉^z} playing the role of the “system fermions.”
The remaining {c^x_ J} and {c^y_ J} operators describe the fermionic modes of the bath, 
while the {b^α_ J} operators introduced above Eq. (<ref>) only participate through their roles in setting the values of the conserved ℤ_2 gauge fields {u_ I,  J}.
The term multiplied by g_ J in Eq. (<ref>) allows fermions to hop between the system and the bath, 
while the term multiplied by B_ J acts on the bath fermions alone. 
Our aim is to prepare the system fermions {c_𝐉^z} in the ground state of the Hamiltonian (<ref>) with g_=0 and no ℤ_2 fluxes, i.e., W_i,j=1 (equivalent to setting u_𝐈,𝐉=1, up to a ℤ_2 gauge transformation).



 §.§ Preparation Protocol

We now present the protocol to prepare the ground state of H_KSL, starting from an arbitrary initial state. Similarly to Refs. <cit.>, the protocol consists of cycles which are repeated until convergence to the ground state is achieved. 
Since the protocol is intended to be implemented on bosonic qubits, in this section we return to the description in terms of the σ and τ spins.

In each cycle, the σ and τ spins are evolved unitarily for a time period T, with a time-dependent Hamiltonian (<ref>) designed to decrease the expectation value of H_KSL. The unitary evolution is followed by a projective measurement of the τ spins in the z basis. The measurement outcomes are used to determine the Hamiltonian of the next cycle. 

In the beginning of the nth cycle, the τ spins are assumed to be in an eigenstate of τ^z_ with known eigenvalues, which we denote by τ_(t_n)=± 1. Before the first cycle, the system is initialized in a state where τ^z_𝐉=+1 for all 𝐉. We then perform the following operations:




    
  * The σ spins are brought to a flux-free state, W_i,j = 1. This can be done by measuring all the flux operators (which are mutually commuting), and annihilating the detected fluxes in pairs by applying appropriate unitary operators <cit.> (see Appendix <ref> for details). 

    
  * The system is evolved from time t_n to t_n+1 = t_n + T with the Hamiltonian H(t-t_n), where H(t) is given in Eq. (<ref>).
    Specifically, we set the effective Zeeman fields to B_𝐉(t) = B(t-t_n)τ_𝐉(t_n)  and the “system-bath couplings” to g_𝐉(t) = g(t-t_n), with the functions g(t), B(t) as shown in Fig. <ref> and given explicitly in Appendix <ref>. 

    
  * The τ^z_ operators are measured, giving new values {τ_(t_n+1)}. A new cycle then begins from step b.  





Importantly, the flux operators W_i,j commute with H(t).

Ideally, they therefore retain their initial values, W_i,j=1. In the presence of decoherence on a real quantum simulator, W_i,j may flip during the cycle. This may require returning to step 1 between cycles, where the flux operators are measured again and unitary operators are applied to correct flipped fluxes, as discussed in Sec. <ref> and Appendix <ref>. 

The idea behind this protocol is that, in analogy to adiabatic demagnetization, the coupling between the σ and τ spins cools the system towards the ground state of H_KSL. 

The τ spins are initially polarized in a direction parallel to the effective Zeeman fields acting on them. 
Adiabatically sweeping the Zeeman fields downward while the coupling g is non-zero tends to transfer energy into the τ spins, decreasing the expectation value of H_KSL. This picture guides the choice of the protocol parameters, B_0, g_1 and T: B_0 should be chosen to be large enough compared to the gap E_gap between the ground state and the lowest excited state of H_KSL, such that sweeping the Zeeman field during the protocol brings the excitations in the system to resonance with the τ spins. The maximum coupling g_1 controls the rate of the cooling. T is chosen such that the time evolution is adiabatic with respect to the system's gap, which sets the maximum allowed values of B_0 and g_1: these are chosen such that max(|B_0|, |g_1|)/T≪ E_gap^2. 

The operation of the cooling protocol is particularly transparent in the fermionic representation, Eq. (<ref>). It is straightforward to check that in the beginning of each cycle, when g_=0, the bath fermions c^x and c^y are decoupled from the system fermions, c^z, and are initialized in the ground state of the bath Zeeman Hamiltonian. Sweeping |B_(t)| downwards induces level crossings between states of the system (c^z) and bath fermions, and when g_ 0, excitations are transferred from the system to the bath. 



 §.§ Performance analysis


The fact that the system and bath can be described as emergent fermions which are coupled to the same gauge field means that a single fermionic excitation can be transferred coherently to the bath. This dramatically accelerates the cooling process compared to simpler “bosonic” adiabatic cooling algorithms of the type described in Refs. <cit.>. In these protocols, only a gauge-neutral pair of fermionic excitations can be transferred from the system to the bath, slowing down the cooling process as the ground state is approached. 

Specifically, the performance of our protocol can be understood from simple considerations. The argument 
mirrors that given in Ref. <cit.> for “gauged cooling” of a one-dimensional quantum Ising model. 
Since fermionic excitations in the system can independently hop into the bath, each excitation has a finite probability to be removed by the bath in each cycle. 
The cooling rate is therefore proportional to the energy density of H_KSL. 

Approximating the cooling process as a continuous time evolution, we obtain a rate equation for the energy density, ε(t):

    ε̇(t) = -C [ε(t)-ε_s],

where C>0 is the cooling rate, and ε_s is the steady state energy density. 
The quantities C and ε_s depend on the protocol parameters, such as T, B_0 and g_1, and in a realistic noisy quantum simulator, also on the noise rate; in the perfectly adiabatic, noiseless case, ε_s→ε_0, the ground state energy density. Solving Eq. (<ref>), we obtain ε(t) = ε_s + [ε(0)-ε_s]e^-C t.
In contrast, in a simple adiabatic cooling protocol where only pairs of fermionic excitations can be transferred to the bath, ε̇∝ -(ε-ε_s)^2, resulting in a power-law approach to the steady state, ε-ε_s ∼ 1/t. 

The fact that the fermionic Hamiltonian (<ref>) is quadratic allows us to analyze the cooling dynamics in a basis in which the excitations are decoupled. Moreover, within the flux-free sector, and assuming that the system is prepared initially such that τ^z_𝐉=1 and hence the signs of the B_𝐉's are spatially uniform, the Hamiltonian can be brought into a translationally invariant form by a unitary transformation. 
This allows an explicit analysis of the performance of the protocol in momentum space, which we present in Appendix <ref>. The analysis verifies the exponential convergence to the steady state with the number of cycles performed, as anticipated above. Moreover, we find that in the ideal (noiseless) case, the steady state energy ε_s depends only on the adiabaticity of the protocol with respect to the system's gap, E_gap, whereas the cooling rate C depends also on the the adiabaticity with respect to the avoided crossing gaps encountered during the time evolution with H(t). In particular, in the limit where T is large compared to 1/E_gap, the ground state can be approached with an exponential accuracy.


While these results are demonstrated for the solvable model H_KSL, which can be mapped to free fermions, we expect them to hold more generally. For example, the performance is not expected to be parametrically affected if quartic interactions between the system fermions are added to Eq. (<ref>). The application of our protocol to general interacting fermionic Hamiltonians is discussed further in Sec. <ref>. 









 

§ NUMERICAL SIMULATIONS

 
 

 §.§ Performance without noise



We simulate the fermion cooling algorithm numerically with and without noise using the efficient procedure described in Appendix <ref> which relies on the model being one of free fermions. In all simulations we choose 𝒥_x=𝒥_y=𝒥_z=𝒥=1, κ=1, in the chiral phase,
 and use protocol parameters B_0=7, g_1=0.5. Periodic boundary conditions are used unless otherwise stated.


 Setting periodic boundary conditions allows us to simulate this translation invariant model in k-space point by point. We cool a system of size 85×85 unit cells in the absence of noise starting from a completely random state. Here, it is assumed that the system has been initialized to the flux-free sector prior to the action of the cooling protocol, which leaves the gauge fields fixed. 
 Time evolution is performed 
 in each momentum sector using an ordinary differential equation solver with a relative tolerance of 10^-6 such that the deviation of the steady state from the ground state is expected to derive predominantly from diabatic transitions.
 
 In the ground state, the Chern number of the c^z_ fermions is equal to 1. To probe the convergence to the ground state, we define the quantity ν^l:

    ν^l=1/2π i∫R^l[∂_k_xR^l,∂_k_yR^l]dk_xdk_y,

 where R^l is the single particle density matrix of either the system (l=sys) or the bath (l=bath) fermions. Specifically, R^sys_s,s'(k)=⟨ c_(k,s)^z† c_(k,s')^z⟩, with an analogous expression for R^bath. 
 Here, c_(k,s)^z is the Fourier transform of c_^z = c_(i, j ,s)^z on sublattice s with wavevector k.  
 In a (pure) Slater determinant state, the quantity ν^l becomes the integrated Chern density. We calculate this value by performing a discrete integral in momentum space.
 
 Fig. <ref> shows ν^sys (solid lines) and ν^bath (dashed lines) at the end of a single cooling cycle before performing the measurement of the τ spins, as a function of the 
 sweep time, T. As the 
 sweep duration is increased, the Chern number of the system approaches 1 and that of the bath approaches -1, even after a single cycle. The fact that our protocol succeeds in preparing a fermionic 
 chiral ground state starting from a 
 
 topologically trivial state
 by a finite depth local unitary circuit is due to 
 
 the chiral state possessing “invertible” topological order; i.e., stacking two systems with opposite Chern numbers yields a trivial state. 
 
 Thus, the topological state of the system's fermions and its inverse (in the bath)
 can be prepared starting from a trivial state [Note, however, that overall, the spin system is in a non-invertible topologically ordered state, since the chiral phase of the KSL has fractional (non-Abelian) statistics.]. 





The inset of the Fig. <ref> shows the energy density after a single cooling cycle as a function of the duration T. As the cycle duration is increased, the system approaches the ground state energy in a single cycle. The deviation from the ground state energy decreases exponentially with increasing T.






 §.§ Performance with noise



In this section, we numerically test the performance of the proposed cooling algorithm in the presence of noise.
For simplicity, the noise is modeled as a uniform depolarizing channel acting on all σ and τ qubits 
at some rate ζ of errors per cooling cycle per qubit. This noise is simulated stochastically as described in Appendix <ref>.

Importantly, 
the fluxes of the gauge field may be excited by the action of noise on the system. In order to overcome this, one should periodically measure and fix these fluxes, while performing our proposed algorithm for cooling the excitations of the fermions {c_ J^z}. Fixing the observed fluxes can be done in an analogous way to the action of error correction in the surface code <cit.> by annihilating them in pairs, as explained in Appendix <ref>. When flux excitations are removed (by applying local spin operations), fermionic excitations may be created, and these are then cooled by subsequent cycles of the cooling protocol. Using the protocol, the energy density of the steady state should be proportional to the noise rate, although there is a finite density of flux excitations in the steady state.

We have simulated a KSL system of size 6× 6 unit cells, cooled in the presence of decoherence for multiple cycles. The results are presented in Fig. <ref>. The energy density is shown at the end of each cycle, with and without the active correction of the flux degrees of freedom.  
In the simulation with flux correction, the fluxes are measured and removed at the end of each cycle. For simplicity, we have assumed that the measurements are ideal, and the flux correction operations are perfect.

In both cases the stochastic action of errors appear as local peaks in the energy density after which the fermionic excitations produced are cooled. However, when errors which excite the flux degrees of freedom occur and these are not corrected, further cycles of fermionic cooling are only able to reach the ground state within the new excited flux sector. After many cycles, a steady state is reached in which the flux degrees of freedom are in a maximally mixed state.





Figure <ref> shows the energy density of the steady state reached by cooling a system of size 10× 10 in the presence of noise while correcting the flux degrees of freedom at each cycle. The energy density is linear in the error rate, demonstrating that cooling of both fermionic and flux excitations is done at a rate proportional to their density.

The active correction of the flux degrees of freedom by measurement and feedback provides a natural tool for further error mitigation by post selection. In cycles which require a large weight Pauli correction in order to annihilate the monitored flux excitations, the error together with its correction likely result in the insertion of fermionic excitations which must be removed by further cooling cycles. We denote the weight of the Pauli correction (defined as the number of gates applied to correct the fluxes) at cycle n by w_n, and calculate its exponential moving average

    w̅_n = ∑_m≤ n w_m e^-n-m/n_0/∑_m≤ n e^-n-m/n_0,

where 1/n_0 is the exponential cooling rate in the absence of noise. We perform post selection in Fig. <ref> by accepting a fraction of the cycles with the smallest values of w̅_n, using n_0=3.
Performing such post selection succeeds in reducing the expectation value of the energy density. In the limit of a very low accepted fraction of the data, the energy density is reduced approximately by a factor of two at low error rates. In this limit, all cycles in which σ spins are affected by errors are removed, but cycles in which τ spins are disturbed remain part of the data. Further post selection can be done by imposing a requirement on the outcomes of the measurements of the τ degrees of freedom from cycle to cycle <cit.>.





§ EXTENSION TO ARBITRARY FERMIONIC HAMILTONIANS




Our protocol for preparing the ground state of the KSL model can be generalized as a method for preparing the ground state of an arbitrary (interacting) fermionic Hamiltonian in 2D. Given some local target Hamiltonian acting on the system Majorana fermions c_𝐉^z, we introduce the gauge fields u_𝐈,𝐉 artificially and modify the target Hamiltonian such that the c_𝐉^z fermions become charged under the gauge field. 
Namely, we modify each hopping term between site 𝐈 and site 𝐉 in the target Hamiltonian by the product of the gauge fields on a path connecting the two sites. 
We then introduce the bath fermions c_𝐉^x,c_𝐉^y, 
and use Eq. (<ref>) to map the fermionic Hamiltonian (including the hopping between the system and bath fermions, used for cooling the system) onto a bosonic local Hamiltonian acting on the σ and τ spins. 

The ground state of the Hamiltonian is prepared according to the protocol above and this state is related to the ground state of the target Hamiltonian by the known gauge transformation as discussed in Sec. <ref>.

Generalizations to other problems, e.g., with complex spinful fermions or multiple orbitals per site, are straightforward. We also note that Hamiltonians with fermion number conservation can be implemented within our protocol; the fermion number changes during each cycle, but the system is cooled towards the ground state that has a well-defined fermion number set by a chemical potential term. 





§ SUMMARY


In this work, we have presented an efficient protocol for preparing certain chiral topological states of matter on quantum simulators. Our prime example is the chiral phase of the Kitaev spin model on the honeycomb lattice. Our protocol can be similarly used to prepare invertible topological phases of fermions, such as Chern insulators and chiral topological superconductors. Importantly, our protocol does not require an a-priori knowledge of the ground state wavefunction, and should perform equally well in the presence of interactions (although the phases we prepare have a description in terms of non-interacting fermions). In the absence of interactions and decoherence, the ground state can be reached after a single cycle, up to adiabatic errors 
(for the non-interacting case, a similar approach was proposed in Ref. <cit.>). 



Our scheme utilizes a simultaneous fermionization of the target system and the cooling bath with which it is coupled, in order to allow the removal of single fermionic excitations from the system to the bath. Crucially, the protocol relies on repeated measurements of the bath spins that detects the fermionic excitations, and local classical feedback that removes them. 

 In practice, for running this algorithm on real quantum hardware, one should choose the time duration T of the unitary evolution and the number of Trotter steps, N_t, to optimize the trade-off between diabatic errors and Trotterization noise (which decrease upon increasing T and N_t) versus the noise deriving from environmental decoherence and low fidelity gates (which are typically proportional to the number of gates applied, and hence to N_t). The optimization can be done even without knowledge of the noise characteristics, by minimizing the measured energy of the steady state with respect to the protocol parameters, such as T and N_t. 



Our protocol offers a parametric advantage compared to other recently-proposed methods to prepare chiral 2D states. Ref. <cit.> describes a method to prepare the chiral state of the Kitaev honeycomb model by adiabatic evolution. Due to the gap closing along the path, the excitation energy of the final state (or the ground state fidelity) scales polynomially with the duration of the time evolution. In contrast, in our protocol, no adiabatic path is required, and the excitation energy scales exponentially with the cycle duration. In Ref. <cit.>, a method for preparing chiral states using sequential quantum circuits is described; this method requires a circuit whose depth scales with the linear size of the system in order to reach a state with a given energy density. Ref. <cit.> provides an entanglement renormalization circuit which can prepare chiral states in a depth that scales logarithmically with system size, however this requires the use of long-range unitary gates. By comparison, within our method the required circuit depth is independent of the system size and local qubit connectivity in 2D is sufficient.



We emphasize that, while our method does not require knowing the ground state in advance, the performance of the protocol depends crucially on the topological phase to which the ground state belongs. Developing practical and efficient methods to prepare ground states in more complicated topological phases, such as fractional Chern insulators, 
is an interesting outstanding challenge. 





G.K. and E.B. were supported by CRC 183 of the Deutsche Forschungsgemeinschaft (subproject A01) and a research grant from
the Estate of Gerald Alexander. M.R. acknowledges the Brown Investigator Award, a program of
the Brown Science Foundation, the University of Washington College of Arts and Sciences, and the Kenneth K.
Young Memorial Professorship for support.





§ REMOVING FLUXES IN 2D


In order to reach the ground state of the 2D KSL, one should prepare the system in a flux-free state of the gauge field u and then proceed to cool the fermionic excitations c^z. 
Fluxes can be removed by measuring the plaquette operators W_i,j which correspond to local 6-body spin operators shown in Eq. (<ref>). 
The excitations W_i,j = -1 can be removed in pairs by products of local unitary operations. We note that the operator σ^α_𝐉 at a vertex of the honeycomb lattice shared between three hexagons anticommutes with two of the three flux operators on these hexagons,

while commuting with the third (as well as all the other plaquettes in the lattice). Therefore, applying such a unitary operator σ^α_𝐉 flips the signs of the two corresponding plaquettes, potentially also inserting fermionic excitations. Choosing some pairing of the measured excited plaquettes, one can construct a string of spin operators ∏_𝐉∈ Sσ^α(𝐉)_𝐉 connecting 
each pair such that the product anti-commutes only with the 
plaquette operators at the end points of the string.


Using the above approach, a finite density of flux exitations can be corrected using measurements and feed-forward unitaries (analogously to the error correction protocol in the surface code <cit.>). The composite error (the original error affecting the fluxes together with the correction) is a product of Pauli operators acting on the σ spins which commutes with all of the flux operators, leaving them unaffected. However, the action of this operator on the system excites a number of residual fermionic excitations proportional to its weight (the number of spins on which it acts). Therefore, assuming a low error rate resulting in a low weight Pauli error with high probability, it is important to construct a correction with a low weight as well by annihilating the flux excitations in local pairs - using a minimum weight perfect matching decoder <cit.>. Otherwise, the density of residual fermionic excitations will scale with system size.

The resulting fermionic excitations left after fusing the fluxes can be removed later by the coupling to the bath. Consequently, in the presence of noise, the combined algorithm (composed of a unitary sweep, measurement of the bath spins and the fluxes, and feed-forward correction of the fluxes) should reach a steady state whose energy density is proportional to the noise rate and is independent of system size. This is the algorithm used in presence of noise, whose results are shown in Figs. <ref> and <ref>.



§ SMOOTH EVOLUTION OF TIME DEPENDENT COUPLINGS



In order to approach the adiabatic limit exponentially with increasing cooling cycle duration T, we choose the time dependence of the couplings B and g to have finite derivatives at all orders. Explicitly, we use a smooth step function 𝒮(x) as used in Ref. <cit.>, defined as

    𝒮(x) = 
    			0,     x<0
    
                1,     x>1
    1/1+e^1/x+1/x-1,    otherwise.

Using this, we express B(t) and g(t) as

    B(t)=B_0[1-𝒮(t/T)],


    g(t)=g_1𝒮(t/t_1)[1-𝒮(t-t_2/T-t_2)].

These curves are shown in Fig. <ref> for t_1=1/4T, t_2=3/4T as chosen in our numerical simulations.



§ ANALYSIS IN MOMENTUM SPACE



In this section we study the cooling dynamics of the quadratic fermionic Hamiltonian (<ref>) in the flux-free sector with uniform τ^z_𝐉=1. We choose a gauge where u_𝐈,𝐉=1, for which the Hamiltonian is translationally invariant,  and Fourier transform the fermionic Hamiltonian to find

    H=∑_kΨ_k^† H_kΨ_k,

where
Ψ_k^†=[ c_(k,A)^z† c_(k,B)^z† c_(k,A)^y† c_(k,B)^y† c_(k,A)^x† c_(k,B)^x† ],









with the 6 × 6 time-dependent Hamiltonian describing the modes at each k given by 







    H_k(t)=([     h_k    G(t)       0;  G^†(t)       0 -iB(t)𝕀;       0  iB(t)𝕀       0 ]).

Here, 𝕀 is a 2× 2 unit matrix, 
h_k is the 2 × 2 submatrix describing the “system” (c^z_(k,s)) fermions, 

    h_k   =[    Δ(k)   if(k); -if(-k)   -Δ(k) ],

and the 2 × 2 matrix G that describes hopping between fermionic system and bath (c^y_(k,s)) is given by

    G(t)=   [ -ig(t)Γ_A         0;         0 -ig(t)Γ_B ],

where we allow the coupling g_ to be sublattice dependent g_(i,j,s)(t) = g(t-t_n)Γ_s rather than spatially uniform as chosen in the main text. 

The functions Δ(k) and f(k) in Eq. (<ref>) are given by

Δ(k)=2κ[sink_x-sink_y+sin(k_y-k_x)]
and f(k) = 𝒥_x+𝒥_ye^-ik_x+𝒥_ye^-ik_y.






 §.§ Adiabatic limit



For the cooling protocol to succeed, it is essential to ensure that for all k the off-diagonal terms G have a non-zero matrix element between the ground state of the target Hamiltonian, h_k, and the ground state of the Hamiltonian Bσ^y of one of the bath modes, c^x_(k,s),c^y_(k,s). Importantly, if the system-bath coupling 
is nonvanishing for all k, 
in the perfectly adiabatic limit  the system reaches its ground state after a single cycle (see discussion in Ref. <cit.>). 
The instantaneous spectrum of H_k in the topological phase with 𝒥_x=𝒥_y=𝒥_z and κ>0 (where the ground state of the system fermions  c^z_(k,s) carries a Chern number of 1 when g=0) is shown at the point k=(2π/3,-2π/3) as a function of B in Fig. <ref>. At this point, we find f(k)=0 and h_k∝σ^z, such that

using only the A sublattice bath (taking Γ_B=0) does not allow cooling of fermionic excitations at this k point.



The existence of this gap closing is due to the non-zero Chern number of the target ground state. As k varies, the spinor corresponding to the system's ground state covers the entire Bloch sphere. Therefore, there must be a value of k for which the system's ground state is orthogonal to the spinor formed by acting with the first column of G on the bath's ground states. At this k point, the coupling to 
the bath on the A sublattice alone cannot cool the system. 


Using two bath sublattices is enough to overcome this problem. As long as matrix G is not singular, the spinors obtained by acting with the columns of G on the system's ground state are not identical, and at least one of these spinors has a non-zero overlap with the system's ground state at each k. Note that, in contrast to the topological phase, in the trivial phase where the Chern number vanishes (obtained, e.g., when 𝒥_z >> 𝒥_x, 𝒥_y and κ=0), one bath sublattice would generically be sufficient (with appropriately chosen coupling). 



 §.§ Diabatic corrections



Deviations from the adiabatic limit have two main effects. First, the protocol
does not converge to its steady state after a single cycle, but rather
converges to the steady state exponentially with the number of cycles
performed. 
Second, the steady state deviates from the ground state of the
system. Importantly, we shall now show that the deviation of the steady
state from the ground state due to diabatic errors is controlled by
the intrinsic energy gap of the system, and not by the avoided crossing
in the middle of the cycle (see Fig. <ref>). The mid-cycle avoided crossing gap controls
the rate of convergence to the steady state.

Returning to the case of Γ_A=Γ_B=1, the diabatic corrections are most straightforward to analyze after
a series of unitary transformations which map the problem at each k into an evolution of three spins.

First, we perform a k-dependent unitary transformation 

    U_k = [ V_k;     V_k;     V_k;     ],

where V_k is chosen such that V_k h_k V_k^† = E_kσ^y,
with E_k=√(|Δ(k)|^2+|f(k)|^2) being the spectrum of the KSL model. 

This brings the Hamiltonian
(<ref>) into a purely antisymmetric form

    H_k=U_kH_kU_k^†=([  E_kσ^y -ig(t)𝕀       0;  ig(t)𝕀       0 -iB(t)𝕀;       0  iB(t)𝕀       0 ]).
















Next, we write the corresponding vector of complex fermionic annihilation operators
Ψ_k=U_kΨ_k in terms of two vectors of Majorana
operators: Ψ_k=1/2(a_k+id_k), where a_k^T=(a_1,k,…,a_6,k)
with a_n,k=a_n,k^†, and similarly for d_k.
Using the antisymmetric property of H_k, the many-body
Hamiltonian decouples when written in terms of a_k, d_k:

    H=1/4∑_k[a_k^TH_ka_k+d_k^TH_kd_k].

For a given k, 
the parts of the Hamiltonian that
depend on a_k and d_k are independent and identical, and it is sufficient to analyze one of them. 





  §.§.§ Mapping to a three spin problem



The part that
depends on a_k 
can be written in terms of three pseudospin operators, 
defined
such that μ_k^z=ia_1,ka_2,k, η_k^z=ia_3,ka_5,k, χ_k^z=ia_4,ka_6,k. 
The basis of the corresponding Hilbert space can be labelled according to the eigenvalues of μ_k^z,η_k^z,χ_k^z; we denote the basis states as |μ_k^z,η_k^z,χ_k^z ⟩. 
The a_k-dependent part of the
Hamiltonian is written as 

    H_k^a(t)=-1/2E_kμ_k^z-1/2B(t)(η_k^z+χ_k^z)+1/2g(t)μ_k^x(η_k^x+χ_k^x).

We have thus mapped the problem into that of three coupled spins. The Hamiltonian of the d_k sector has
the same form, and can be treated similarly. 
Note that the μ_k,η_k,χ_k spins are distinct from the physical ones, σ_J and τ_J.



  §.§.§ Solution of the three spin problem



At g=0 the Hamiltonian is diagonalized in the orthonormal basis 1/√(2)(|± 1,1,-1⟩-|± 1,-1,1⟩), |± 1,1,1⟩, 1/√(2)(|± 1,1,-1⟩+|± 1,-1,1⟩), |± 1,-1,-1⟩ with energies given by ∓ E_k/2, ∓ E_k/2-B, ∓ E_k/2, ∓ E_k/2+B. At finite g, the first two states remain decoupled. The Hamiltonian in the complementary subspace written in this basis is given by


    H=[ -E_k/2-B                      √(2)g         ;           E_k/2-B    √(2)g                           ;             √(2)g   -E_k/2                      √(2)g;    √(2)g                      E_k/2    √(2)g         ;                               √(2)g -E_k/2+B         ;                      √(2)g                    E_k/2+B;          ],

and its spectrum is shown as a function of B for g>0 (solid lines) and for g=0 (dotted lines) in Fig. <ref>b. Turning on a non-zero g couples 1/√(2)(|± 1,1,-1⟩+|± 1,-1,1⟩) to |∓ 1,1,1⟩ and |∓ 1,-1,-1⟩. 
(Note that [μ_k^z η_k^z χ_k^z, H^a_k]=0: the even and odd parity sectors of total z pseudospin are decoupled.) 



Each cooling cycle is initialized with the bath spins η_k, χ_k in their respective ground states, i.e., the system is initialized in some mixture of the states |± 1,1,1⟩. For g=0, the energies of the eigenstates |-1,1,1⟩ and 1/√(2)(| 1,1,-1⟩+| 1,-1,1⟩) cross at B=E_k, while the state |1,1,1⟩ remains spectrally separated from the only state it is directly coupled to, 1/√(2)(| -1,1,-1⟩+| -1,-1,1⟩). With a finite g, the crossing becomes an avoided crossing and in the adiabatic limit, the states |-1,1,1⟩ and 1/√(2)(| 1,1,-1⟩+| 1,-1,1⟩) are exchanged under the time evolution, while | 1,1,1 ⟩ 
evolves trivially.

Away from the adiabatic limit, we denote the diabatic transition probability for the | -1, 1, 1⟩, 1/√(2)(| 1,1,-1⟩+| 1,-1,1⟩) pair as P_1,k and the diabatic transition probability 
between | 1,1,1⟩ and 1/√(2)(| -1,1,-1⟩+| -1,-1,1⟩) as P_2,k (see Fig. <ref>b).
For g(t)/E_k≪ 1 and assuming g(t) is constant for times where B(t) ≈ E_k, P_1,k is approximately given by the Landau-Zener formula, P_1,k≈ e^-c_1 Tg^2/B_0 for some constant c_1. The second transition probability, 
P_2,k, is determined by the single particle gap of the target Hamiltonian, P_2,k≈ e^-c_2 TE_k^2/B_0≥ e^-c_2 TE_gap^2/B_0, where E_gap is the minimal gap. We neglect the possibility for transitions to the states |± 1,-1,-1⟩ as a second order effect, occurring with probability of the order of P_1,kP_2,k.



  §.§.§ Analysis of the cooling dynamics






































 A single cycle of the cooling protocol, including the measurement of the ancilla spins, can be described as a map from the initial density matrix for the system (corresponding to the pseudospin μ_k) to the final one, i.e., the dynamics realize a quantum channel. 
The protocol consists of applying this quantum channel to the system multiple times, until a steady state is reached.

We number the four states participating in the dynamics | 1,1,1⟩, 1/√(2)(| 1,1,-1⟩+| 1,-1,1⟩), | -1, 1, 1⟩, 1/√(2)(| -1,1,-1⟩+| -1,-1,1⟩) by 0,1,2,3 for brevity and denote the corresponding density matrix by

    ρ=[ ρ_00 ρ_01 ρ_02 ρ_03; ρ_10 ρ_11 ρ_12 ρ_13; ρ_20 ρ_21 ρ_22 ρ_23; ρ_30 ρ_31 ρ_32 ρ_33 ].


The final step of the protocol, in which τ_J^z are measured and the signs of B_J are adjusted accordingly for the next cycle, is equivalent to a resetting operation taking η_k^z,χ_k^z→1 after the measurement. The action of this operation on the density matrix is given by

    ρreset→[ ρ_00+ρ_11         0 ρ_02+ρ_13         0;         0         0         0         0; ρ_20+ρ_31         0 ρ_22+ρ_33         0;         0         0         0         0 ].




The unitary time evolution step acts as

    ρunitary→𝒰ρ𝒰^†,

where

    𝒰=([            √(P_2,k)                                              e^iθ_1√(P_2,k);                               -√(P_1,k)      e^iθ_2√(P_1,k);                          e^iϕ_2√(P_1,k) e^iθ_2+iϕ_2√(P_1,k);     -e^iϕ_1√(P_2,k)                                         e^iθ_1+iϕ_1√(P_2,k) ]),

where P_1,k and P_2,k are the transition probabilities between the different states (see Sec. <ref> and Fig. <ref>b), and
we denote the complementary probabilities to P_n,k by P_n,k=1-P_n,k. Finally, under the full protocol cycle consisting of unitary evolution followed by measurement and reset,

    [ ρ'_00; ρ'_02; ρ'_20; ρ'_22 ]=M[ ρ_00; ρ_02; ρ_20; ρ_22 ],

where

    M=[          P_2,k                                        P_1,k;                 √(P_2,kP_1,k) -√(P_2,kP_1,k);                -√(P_2,kP_1,k)  √(P_1,kP_2,k);          P_2,k                                        P_1,k ].

Here, the complex phases
e^iθ_1,e^iθ_2,e^iϕ_1,e^iϕ_2 cancel. If we take P_2,k→0, then the steady state of the protocol is the desired ground state (μ_k^z=1) since there is no probability
for any transition.
In general, the steady state of these dynamics is given by 

    M[ P_1,k/P_1,k+P_2,k;                 0;                 0; P_2,k/P_1,k+P_2,k ]=[ P_1,k/P_1,k+P_2,k;                 0;                 0; P_2,k/P_1,k+P_2,k ];


the steady state is a mixture of the ground state (μ_k^z=1) and the excited state (μ_k^z=-1). The probability to be in the ground state is given by

    P_ steady, k =  1-P_1,k/1-P_1,k+P_2,k.


We denote the eigenvalues of the quantum channel by λ_n,k (where λ_1,k=1 corresponds to the steady state). 
The eigenvalues, λ_2,k, λ_3,k, and λ_4,k 
determine the rate in cycles at which the steady state is exponentially approached. These are given by

    λ_2,k=√(1-P_2,k)√(P_1,k)+√(P_2,k)√(1-P_1,k),


    λ_3,k=√(1-P_2,k)√(P_1,k)-√(P_2,k)√(1-P_1,k),


    λ_4,k=P_1,k-P_2,k.

The largest among these (i.e., the eigenvalue closest to 1) is λ_2,k, and hence this eigenvalue controls the approach to the steady state: asymptotically, the deviation from the steady state scales with the number of cycles n as e^-n/n_0, where n_0 = 1/ln(1/λ_2,k)

[Note that, in general, λ_2,k≤ 1. In the special case P_1,k=1-P_2,k, λ_2,k=1 and the steady state becomes
degenerate. For any other value of P_1,k, the system converges to the steady state.]. 








Interestingly, Eq. (<ref>) shows that for P_2,k→0, the steady state is exactly the ground state, 
independent of the value of P_1,k.
This is due to the fact that, for P_2,k=0, there is no way to escape from the ground state once it is reached. However, the rate of approach to the ground state does depend on P_1,k. 
















































§ SINGLE-PARTICLE DENSITY MATRIX SIMULATION OF FREE FERMIONS



The free evolution of a quantum many-body system under a quadratic fermionic Hamiltonian, including measurements and resets of single sites and Pauli errors, can be traced numerically efficiently, requiring only 
time and memory that scale as polynomials of the system size. Here, we briefly describe the technique used to simulate the noisy dynamics of the cooling process discussed in Sec. <ref>. The technique is the same as that used in Ref. <cit.>.

We keep track of the single particle density matrix R_𝐈,𝐉^α,β=⟨ ic_𝐈^αc_𝐉^β⟩ as well as the values of the ℤ_2 gauge field u_𝐈,𝐉=± 1 as they evolve under the noisy cooling dynamics.





This allows us to calculate any quadratic observables, and the energy in particular.







During free evolution under the Hamiltonian (<ref>) of the main text within a fixed sector of the gauge field u, the single particle density matrix evolves under a unitary matrix W(t) by conjugation, R(t)=W^†(t)R(0)W(t). This unitary is given by the time ordered exponentiation of the Hamiltonian in matrix form, and can be found numerically by Trotterization or using an ordinary differential equation solver. In Sec. <ref> we use second order Trotterization, making N_t steps per cooling cycle. Importantly, the unitary depends on the values of the gauge field which are looked up for its computation.

The resetting operation of the bath fermions at the end of each cycle ic_𝐉^yc_𝐉^x→ 1 is handled as in Ref. <cit.>. This amounts to setting all matrix elements of R_𝐈,𝐉^α,β in the two rows and the two columns corresponding to 𝐉,x and to 𝐉,y to zero, and then setting the 2×2 block at their intersection to the values found at ic_𝐉^yc_𝐉^x=1.


























































































Depolarizing noise is implemented through stochastic application of Pauli errors before each Trotter step of the time evolution and for each spin ({σ_𝐉,τ_𝐉}) with probability of ζ/N_t. The Pauli τ_𝐉^α operators are bilinear in the c fermions (<ref>), such that they correspond to evolution under a quadratic fermionic Hamiltonian. Pauli errors σ_𝐉^α act on the b Majorana fermions, flipping the signs of two of the three gauge fields incident to site 𝐉. 

Measuring the flux degrees of freedom within this framework simply requires reading off the product of the values of the ℤ_2 gauge fields around the corresponding plaquettes. These values are well defined at every point in time for each trajectory in these simulations, and their fluctuations can be obtained by averaging over multiple trajectories with stochastic noise. Correcting these fluxes in order to return to the flux-free sector, as described in Appendix <ref>, is done by applying a unitary correction given by a product of Pauli σ_𝐉^α operators. 


















