Efficiently preparing chiral states via fermionic cooling on bosonic quantum hardware

Simulating quantum many-body systems, both in and out of equilibrium, is one of most promising applications of near-term quantum computers^1,2. 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³. These unitary circuit algroithms include variational methods^4,5,6, adiabatic processes^7,8,9,10, and effective imaginary time evolution^11,12,13. In addition, protocols for cooling the system by coupling it to a simulated thermal bath have been proposed^{14,15,16,17,18,19,20,21}. 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^22,23. The principle is illustrated by the example of the toric code model²⁴: 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²⁵. The KSL model can be mapped into a system of emergent Majorana fermions coupled to a static ({{mathbb{Z}}}_{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²⁶. 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 ({{mathbb{Z}}}_{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. ²⁷ to two dimensions. (By bosonic qubits, we mean that operators belonging to different qubits commute with each other. Qubits can be thought of as hard-core bosons). The bath fermions are then used to extract energy and entropy from the system, by performing simulated cyclic adiabatic cooling^20,27 followed by measurements of the auxiliary qubits. The steps of the protocol are illustrated in Fig. 1.

The system’s fermionic degrees of freedom (black) together with a fermionic “bath” (white) are encoded using bosonic qubits such that the two are charged under the same ({{mathbb{Z}}}_{2}) gauge field. a The system is brought to a flux-free state by measuring the flux operators and annihilating defects in pairs by applying a string of spin operators. b The system is evolved unitarily with the time-dependent Hamiltonian in Eq. (2), allowing single fermionic excitations to hop coherently from the system to the bath. c At the end of the cooling cycle, the hopping between the system and bath sites is turned off. The location of the fermionic excitations in the bath is measured and the measurement outcomes are used to determine the Hamiltonian for the next cycle.

Full size image

Our protocol prepares the ground state of the chiral phase parametrically more efficiently than other methods, such as adiabatic preparation⁷ or ‘naive’ adiabatic cooling using a simulated bosonic bath^15,20. 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²⁰).

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 ({{mathbb{Z}}}_{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.

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-(frac{1}{2}) degrees of freedom on a honeycomb lattice with strongly anisotropic exchange interactions^25,28,29. The Hamiltonian is written as:

Here, J = (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 ({sigma }_{{{bf{J}}}}^{alpha }) denotes the α = x, y, z Pauli matrix at site J. The bonds of the lattice are partitioned into three sets – x, y, and z bonds –according to their geometric orientations (see Fig. 2). We denote the corresponding Kitaev anisotropic exchange couplings by ({{{mathcal{J}}}}_{x},{{{mathcal{J}}}}_{y},) and ({{{mathcal{J}}}}_{z}), respectively. The κ term is a three-spin interaction acting on three neighboring sites denoted by (leftlangle right.leftlangle right.{{bf{J}}},{{bf{K}}},{{bf{L}}}left.rightrangle left.rightrangle), such that both J and K are nearest neighbors of site L. The κ term breaks time reversal symmetry, and opens a gap in the bulk that drives the system into the chiral non-Abelian phase²⁵.

Each unit cell consists of two system spins (black) and two bath spins (white). The bonds are colored according to orientation: x-blue, y-green, z-red. The five-spin term κ drives the system into the chiral phase. A product of the unitaries ({{{mathcal{U}}}}_{{{bf{I}}},{{bf{J}}}}) along a path removes a pair of fermionic excitations from the bath. Instead of applying this non-local operation to the state, any subsequent operator acting on the system is conjugated by it.

Full size image

In order to facilitate an efficient protocol to prepare the ground state of (1), we introduce a modified KSL Hamiltonian, ({widetilde{H}}_{{{rm{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. 2, and is composed of a sum of two terms:

where ({widetilde{H}}_{{{rm{KSL}}}}) is the modified KSL Hamiltonian:

The ({tau }_{{{bf{J}}}}^{z}) operators all commute with ({widetilde{H}}_{{{rm{KSL}}}}). Clearly, ({widetilde{H}}_{{{rm{KSL}}}}) is identical to H_KSL in Eq. (1) in the sector ({tau }_{{{bf{J}}}}^{z}=+1). In fact, the two Hamiltonians are related by a unitary transformation for any state of the τ spins with even parity, ({prod }_{{{bf{J}}}}{tau }_{{{bf{J}}}}^{z}=1). To see this, we note that all of the unitary operators ({{{mathcal{U}}}}_{{{bf{I}}},{{bf{J}}}}={tau }_{{{bf{I}}}}^{x}{sigma }_{{{bf{I}}}}^{alpha }{sigma }_{{{bf{J}}}}^{alpha }{tau }_{{{bf{J}}}}^{x}), for 〈I, J〉 ∈ α-bonds, commute with ({widetilde{H}}_{{{rm{KSL}}}}) and anticommute individually with ({tau }_{{{bf{I}}}}^{z}) and ({tau }_{{{bf{J}}}}^{z}). The transformation from a sector with a given ({{tau }_{{{bf{J}}}}^{z}}) of even parity to the sector {({tau }_{{{bf{J}}}}^{z}=1)} can be constructed by annihilating pairs of ({tau }_{{{bf{K}}}}^{z},{tau }_{{{bf{L}}}}^{z}=-1) using a product of the unitaries ({{{mathcal{U}}}}_{{{bf{I}}},{{bf{J}}}}) along a path connecting the points K, L. Thus, preparing the ground state of ({widetilde{H}}_{{{rm{KSL}}}}) in any known even-parity sector of ({{tau }_{{{bf{J}}}}^{z}}) is equivalent to preparing the ground state of the original KSL Hamiltonian, Eq. (1).

The “control” Hamiltonian H_c(t), which is used to cool into the ground state of ({widetilde{H}}_{{{rm{KSL}}}}), is chosen as a time-dependent effective Zeeman field acting on τ_J with (hat{x}) and (hat{z}) components g_J(t) and B_J(t), respectively,

The time dependence of B_J(t) and g_J(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}_{{{bf{J}}}}^{x},{c}_{{{bf{J}}}}^{y},{c}_{{{bf{J}}}}^{z},{b}_{{{bf{J}}}}^{x},{b}_{{{bf{J}}}}^{y},{b}_{{{bf{J}}}}^{z}}) per site, subject to the constraint

The spin operators are related to the fermionic ones by:

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²⁵. 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. (5) and (6), we can write

Using this mapping, we express Hamiltonian (2) as

where we identify the set of conserved ({{mathbb{Z}}}_{2}) gauge fields

The gauge-invariant flux W_i,j of the ({{mathbb{Z}}}_{2}) gauge field is defined on each hexagonal plaquette that corresponds to the unit cell (i, j), as the product of u_I,J on its edges. In terms of the spin degrees of freedom, this can be written as

The fermionic Hamiltonian (8) with g_J = 0, B_J = 0 is identical to the Hamiltonian found by fermionizing the original KSL model²⁵, with ({{c}_{{{bf{J}}}}^{z}}) playing the role of the “system fermions.” The remaining ({{c}_{{{bf{J}}}}^{x}}) and ({{c}_{{{bf{J}}}}^{y}}) operators describe the fermionic modes of the bath, while the ({{b}_{{{bf{J}}}}^{alpha }}) operators introduced above Eq. (5) only participate through their roles in setting the values of the conserved ({{mathbb{Z}}}_{2}) gauge fields {u_I,J}. The term multiplied by g_J in Eq. (8) 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}_{{{bf{J}}}}^{z}}) in the ground state of the Hamiltonian (8) with g_J = 0 and no ({{mathbb{Z}}}_{2}) fluxes, i.e., W_i,j = 1 (equivalent to setting u_I,J = 1, up to a ({{mathbb{Z}}}_{2}) gauge transformation).

Preparation Protocol

We now present the protocol to prepare the ground state of ({widetilde{H}}_{{{rm{KSL}}}}), starting from an arbitrary initial state. Similarly to Refs. ^20,27, 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 (2) designed to decrease the expectation value of ({widetilde{H}}_{{{rm{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 ({tau }_{{{bf{J}}}}^{z}) with known eigenvalues, which we denote by τ_J(t_n) = ± 1. Before the first cycle, the system is initialized in a state where ({tau }_{{{bf{J}}}}^{z}=+1) for all J. We then perform the following operations:

1. a.

   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³⁰ (see Methods subsection “Removing fluxes in 2D” for details).

2. b.

   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. (2). Specifically, we set the effective Zeeman fields to B_J(t) = B(t − t_n)τ_J(t_n) and the “system-bath couplings” to g_J(t) = g(t − t_n), with the functions g(t), B(t) as shown in Fig. 3 and given explicitly in Methods subsection “Smooth evolution of time dependent couplings”.

   

   The time dependence of the couplings g and B of the control Hamiltonian [Eq. (4)]. Explicit expressions are given in Methods subsection “Smooth evolution of time dependent couplings”.

   Full size image

3. c.

   The ({tau }_{{{bf{J}}}}^{z}) operators are measured, giving new values {τ_J(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. “Performance with noise” and in Methods subsection “Removing fluxes in 2D”.

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 ({widetilde{H}}_{{{rm{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 ({widetilde{H}}_{{{rm{KSL}}}}). This picture guides the choice of the protocol parameters, B₀, g₁ and T: B₀ should be chosen to be large enough compared to the gap E_gap between the ground state and the lowest excited state of ({widetilde{H}}_{{{rm{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₁ 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₀ and g₁: these are chosen such that (,{rm{max}}(| {B}_{0}| ,| {g}_{1}| )/Tll {E}_{{rm{gap}},}^{2}).

The operation of the cooling protocol is particularly transparent in the fermionic representation, Eq. (8). It is straightforward to check that in the beginning of each cycle, when g_J = 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_J(t)∣ downwards induces level crossings between states of the system (c^z) and bath fermions, and when g_J ≠ 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. ^{14,15,16,17,18,19,20,21}. 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. ²⁷ 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 ({widetilde{H}}_{{{rm{KSL}}}}). Approximating the cooling process as a continuous time evolution, we obtain a rate equation for the energy density, ε(t):

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₀ and g₁, and in a realistic noisy quantum simulator, also on the noise rate; in the perfectly adiabatic, noiseless case, ε_s → ε₀, the ground state energy density. Solving Eq. (11), we obtain ε(t) = ε_s + [ε(0) − ε_s]e^−Ct. In contrast, in a simple adiabatic cooling protocol where only pairs of fermionic excitations can be transferred to the bath, (dot{varepsilon }propto -{(varepsilon -{varepsilon }_{{{rm{s}}}})}^{2}), resulting in a power-law approach to the steady state, ε − ε_s ~ 1/t.

The fact that the fermionic Hamiltonian (8) 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 ({tau }_{{{bf{J}}}}^{z}=1) and hence the signs of the B_J’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 Supplementary Note 1, Supplementary Fig. 1 and Supplementary Fig. 2. 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 ({widetilde{H}}_{{{rm{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. (3). In the presence of interactions, excitations retain a finite overlap with single fermions since the non-interacting model is gapped. The application of our protocol to general interacting fermionic Hamiltonians is discussed further in Sec. “Extension to arbitrary fermionic Hamiltonians”.

Numerical simulations

Performance without noise

We simulate the fermion cooling algorithm numerically with and without noise using the efficient procedure described in Methods subsection “Single-particle density matrix simulation of free fermions” which relies on the model being one of free fermions. In all simulations we choose ({{{mathcal{J}}}}_{x}={{{mathcal{J}}}}_{y}={{{mathcal{J}}}}_{z}={{mathcal{J}}}=1), κ = 1, in the chiral phase, and use protocol parameters B₀ = 7, g₁ = 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⁻⁶ 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}_{{{bf{J}}}}^{z}) fermions is equal to 1. To probe the convergence to the ground state, we define the quantity ν^l:

where R^l is the single particle density matrix of either the system (l = sys) or the bath (l = bath) fermions. Specifically, ({R}_{s,s{prime} }^{{{rm{sys}}}}({{boldsymbol{k}}})=langle {c}_{({{boldsymbol{k}}},s)}^{z{dagger} }{c}_{({{boldsymbol{k}}},s^{prime} )}^{z}rangle), with an analogous expression for R^bath. Here, ({c}_{({{boldsymbol{k}}},s)}^{z}) is the Fourier transform of ({c}_{{{bf{J}}}}^{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.

Figure 4 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.

A KSL system of size 85 × 85 with ({{mathcal{J}}}=1,kappa =1) starting from a randomly initialized state in the flux-free sector is cooled for a single cycle using protocol parameters g₁ = 0.5, B₀ = 7. a A finite size version of the spectral Chern number (12) of the system (solid lines) and the bath (dashed lines) before performing the reset is shown as a function of the cycle duration T. The inset shows the energy density vs. T. In the adiabatic limit, the system approaches the ground state, the Chern number of the system approaches 1 and that of the bath approaches  − 1 in a single cycle as illustrated in b.

Full size image

The inset of the Fig. 4 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 Methods subsection “Single-particle density matrix simulation of free fermions”.

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}_{{{bf{J}}}}^{z}}). Fixing the observed fluxes can be done in an analogous way to the action of error correction in the surface code²⁶ by annihilating them in pairs, as explained in Methods subsection “Removing fluxes in 2D”. 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. 5. 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.

A KSL system of size 6 × 6 unit cells initialized in its ground state is cooled in the presence of a decoherence rate of ζ = 10⁻³ errors per cycle per qubit. The cycle duration is set to T = 20 and N_t = 800 Trotter steps are used per cycle. The energy density is plotted at the end of each cycle. The case where fluxes are corrected every cycle is shown in blue and the case where fluxes are left uncorrected is in orange.

Full size image

Figure 6 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 energy density of a KSL model of size 10 × 10 unit cells in the steady state reached after many cooling cycles as a function of the rate of depolarizing noise. The energies are averaged over 2000 cycles per data point and error bars represent the standard error of the mean determined for an effective sample size of 2000/(2n₀) taking into account exponentially decaying correlations over approximately n₀ = 3 cycles. The cycle duration is set to T = 10 and N_t = 400 Trotter steps are used per cycle. The flux degrees of freedom are monitored and corrected at the end of each cycle. Post selection is performed with respect to the weight of decoded Pauli correction applied to fix the detected flux excitations, exponentially averaged over recent cycles ({bar{w}}_{n}). Cycles requiring Pauli corrections of high weight are omitted. Different colors correspond to different fractions of cycles accepted as ranked by this criterion. The data is fit by linear functions (dashed lines).

Full size image

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

where 1/n₀ is the exponential cooling rate in the absence of noise. We perform post selection in Fig. 6 by accepting a fraction of the cycles with the smallest values of ({bar{w}}_{n}), using n₀ = 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²⁰.

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}_{{{bf{J}}}}^{z}), we introduce the gauge fields u_I,J artificially and modify the target Hamiltonian such that the ({c}_{{{bf{J}}}}^{z}) fermions become charged under the gauge field. Namely, we modify each hopping term between site I and site J 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}_{{{bf{J}}}}^{x},{c}_{{{bf{J}}}}^{y}), and use Eq. (6) 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. “Preparing a chiral spin liquid”.

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.

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. ²³).

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. ³⁰ 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. ³¹, 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. ²³ proposes to prepare a “quantum double” state composed of the chiral state of interest and its time-reversed partner. The doubled state supports a gapped boundary, and can be prepared using measurements and feedback in one shot. Then, a finite-depth quasi-local unitary evolution can decouple the state from its time reversal partner. This idea has some similarities to the method we propose both in terms of the use of non-unitary operations and the scaling achieved. However, the routes are distinct and, in particular, here we give an explicit detailed construction of the protocol whereas in Ref. ²³ the idea is presented and justified on general theoretical grounds without providing an explicit protocol for the unitary evolution step (as would be needed to enable an implementation in practice).

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.

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. (10). The excitations W_i,j = − 1 can be removed in pairs by products of local unitary operations. We note that the operator ({sigma }_{{{bf{J}}}}^{alpha }) 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 ({sigma }_{{{bf{J}}}}^{alpha }) 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 ({prod }_{{{bf{J}}}in S}{sigma }_{{{bf{J}}}}^{alpha ({{bf{J}}})}) 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²⁶). 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³². 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. 5 and 6.

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 ({{mathcal{S}}}(x)) as used in Ref. ²⁷, defined as

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

These curves are shown in Fig. 3 for ({t}_{1}=frac{1}{4}T), ({t}_{2}=frac{3}{4}T) as chosen in our numerical simulations.

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. “Performance with noise”. The technique is the same as that used in Ref. ²⁷.

We keep track of the single particle density matrix ({R}_{{{bf{I}}},{{bf{J}}}}^{alpha ,beta }=langle i{c}_{{{bf{I}}}}^{alpha }{c}_{{{bf{J}}}}^{beta }rangle) as well as the values of the ({{mathbb{Z}}}_{2}) gauge field u_I,J = ±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 Eq. (8) 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. “Performance with noise” of the main text 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 (i{c}_{{{bf{J}}}}^{y}{c}_{{{bf{J}}}}^{x}to 1) is handled as in Ref. ²⁷. This amounts to setting all matrix elements of ({R}_{{{bf{I}}},{{bf{J}}}}^{alpha ,beta }) in the two rows and the two columns corresponding to J, x and to J, y to zero, and then setting the 2 × 2 block at their intersection to the values found at (i{c}_{{{bf{J}}}}^{y}{c}_{{{bf{J}}}}^{x}=1).

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

Measuring the flux degrees of freedom within this framework simply requires reading off the product of the values of the ({{mathbb{Z}}}_{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 Methods subsection “Removing fluxes in 2D”, is done by applying a unitary correction given by a product of Pauli ({sigma }_{{{bf{J}}}}^{alpha }) operators.