Efficiently preparing chiral states via fermionic cooling on bosonic quantum hardware

Introduction

Simulating quantum many-body systems, both in and out of equilibrium, is one of most promising applications of near-term quantum computers1,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 circuit3. These unitary circuit algroithms include variational methods4,5,6, adiabatic processes7,8,9,10, and effective imaginary time evolution11,12,13. In addition, protocols for cooling the system by coupling it to a simulated thermal bath have been proposed14,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 feedback22,23. The principle is illustrated by the example of the toric code model24: 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 excitations25. 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 detected26. 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. 27 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 cooling20,27 followed by measurements of the auxiliary qubits. The steps of the protocol are illustrated in Fig. 1.

Fig. 1: Scheme for simulated cooling of a fermionic system using bosonic qubits.
figure 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 preparation7 or ‘naive’ adiabatic cooling using a simulated bosonic bath15,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 rate20).

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.

Results

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 interactions25,28,29. The Hamiltonian is written as:

$${H}_{{{rm{KSL}}}} = -sum limits_{alpha in {x,y,z}}{{{mathcal{J}}}}_{alpha }sum limits_{leftlangle right.{{bf{I}}},{{bf{J}}}left.rightrangle in alpha ,{rm{-bonds}},}{sigma }_{{{bf{I}}}}^{alpha }{sigma }_{{{bf{J}}}}^{alpha }\ -kappa sumlimits_{leftlangle right.leftlangle right.{{bf{J}}},{{bf{K}}},{{bf{L}}}left.rightrangle left.rightrangle }{sigma }_{{{bf{J}}}}^{x}{sigma }_{{{bf{K}}}}^{y}{sigma }_{{{bf{L}}}}^{z}.$$
(1)

Here, J = (ijs) denotes the sites of the honeycomb lattice, where i {1, …, Nx} and j {1, …, Ny} label the unit cell and s {AB} labels the two sublattices, and ({sigma }_{{{bf{J}}}}^{alpha }) denotes the α = xyz 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 phase25.

Fig. 2: The generalized KSL model.
figure 2

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:

$$H(t)={widetilde{H}}_{{{rm{KSL}}}}+{H}_{{{rm{c}}}}(t),$$
(2)

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

$${widetilde{H}}_{{{rm{KSL}}}}= -mathop{sum}limits_{alpha in {x,y,z}}{{{mathcal{J}}}}_{alpha } mathop{sum}limits_{langle {{bf{I}}},{{bf{J}}}rangle in alpha ,{rm{-bonds}},}{sigma }_{{{bf{I}}}}^{alpha }{sigma }_{{{bf{J}}}}^{alpha }{tau }_{{{bf{I}}}}^{z}{tau }_{{{bf{J}}}}^{z}\ -kappa mathop{sum}limits_{langle langle {{bf{J}}},{{bf{K}}},{{bf{L}}}rangle rangle }{sigma }_{{{bf{J}}}}^{x}{tau }_{{{bf{J}}}}^{z}{sigma }_{{{bf{K}}}}^{y}{tau }_{{{bf{K}}}}^{z}{sigma }_{{{bf{L}}}}^{z}.$$
(3)

The ({tau }_{{{bf{J}}}}^{z}) operators all commute with ({widetilde{H}}_{{{rm{KSL}}}}). Clearly, ({widetilde{H}}_{{{rm{KSL}}}}) is identical to HKSL 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 〈IJ α-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 KL. 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 Hc(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 gJ(t) and BJ(t), respectively,

$${H}_{{{rm{c}}}}(t)=-mathop{sum}limits_{{{bf{J}}}}{B}_{{{bf{J}}}}(t){tau }_{{{bf{J}}}}^{z}-mathop{sum}limits_{{{bf{J}}}}{g}_{{{bf{J}}}}(t){tau }_{{{bf{J}}}}^{x}.$$
(4)

The time dependence of BJ(t) and gJ(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

$${c}_{{{bf{J}}}}^{x}{c}_{{{bf{J}}}}^{y}{c}_{{{bf{J}}}}^{z}{b}_{{{bf{J}}}}^{x}{b}_{{{bf{J}}}}^{y}{b}_{{{bf{J}}}}^{z}=i.$$
(5)

The spin operators are related to the fermionic ones by:

$${sigma }_{{{bf{J}}}}^{alpha }=-frac{i}{2}mathop{sum}limits_{beta ,gamma }{varepsilon }^{alpha beta gamma }{b}_{{{bf{J}}}}^{beta }{b}_{{{bf{J}}}}^{gamma },quad {tau }_{{{bf{J}}}}^{alpha }=-frac{i}{2} mathop{sum}limits_{beta ,gamma }{varepsilon }^{alpha beta gamma }{c}_{{{bf{J}}}}^{beta }{c}_{{{bf{J}}}}^{gamma },$$
(6)

where αβγ {xyz} and εαβγ is the totally antisymmetric tensor. This transformation of the spins to fermions is analogous to the one used by Kitaev25. 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

$${tau }_{{{bf{J}}}}^{z}{sigma }_{{{bf{J}}}}^{alpha }=-i{c}_{{{bf{J}}}}^{z}{b}_{{{bf{J}}}}^{alpha }.$$
(7)

Using this mapping, we express Hamiltonian (2) as

$$begin{array}{rcl}H(t)&=&sumlimits_{alpha in {x,y,z}}{{{mathcal{J}}}}_{alpha }sumlimits_{leftlangle right.{{bf{I}}},{{bf{J}}}left.rightrangle in alpha ,{rm{-bonds}},}{u}_{{{bf{I}}},{{bf{J}}}}i{c}_{{{bf{I}}}}^{z}{c}_{{{bf{J}}}}^{z}\ &&-ikappa sumlimits_{leftlangle right.leftlangle right.{{bf{J}}},{{bf{K}}},{{bf{L}}}left.rightrangle left.rightrangle }{u}_{{{bf{J}}},{{bf{L}}}}{u}_{{{bf{L}}},{{bf{K}}}}{c}_{{{bf{J}}}}^{z}{c}_{{{bf{K}}}}^{z}\ &&-sumlimits_{{{bf{J}}}}{B}_{{{bf{J}}}}(t)i{c}_{{{bf{J}}}}^{y}{c}_{{{bf{J}}}}^{x}-sumlimits_{{{bf{J}}}}{g}_{{{bf{J}}}}(t)i{c}_{{{bf{J}}}}^{z}{c}_{{{bf{J}}}}^{y},end{array}$$
(8)

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

$$left{{u}_{{{bf{I}}},{{bf{J}}}}=i{b}_{{{bf{I}}}}^{alpha }{b}_{{{bf{J}}}}^{alpha }| leftlangle right.{{bf{I}}},{{bf{J}}}left.rightrangle in alpha ,{rm{-bonds}},,alpha in left{x,y,zright}right}.$$
(9)

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

$${W}_{i,j}={sigma }_{(i,j,B)}^{x}{sigma }_{(i+1,j,A)}^{z}{sigma }_{(i+1,j,B)}^{y}{sigma }_{(i+1,j+1,A)}^{x}{sigma }_{(i,j+1,B)}^{z}{sigma }_{(i,j+1,A)}^{y}.$$
(10)

The fermionic Hamiltonian (8) with gJ = 0, BJ = 0 is identical to the Hamiltonian found by fermionizing the original KSL model25, 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 {uI,J}. The term multiplied by gJ in Eq. (8) allows fermions to hop between the system and the bath, while the term multiplied by BJ 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 gJ = 0 and no ({{mathbb{Z}}}_{2}) fluxes, i.e., Wi,j = 1 (equivalent to setting uI,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(tn) = ± 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, Wi,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 operators30 (see Methods subsection “Removing fluxes in 2D” for details).

  2. b.

    The system is evolved from time tn to tn+1 = tn + T with the Hamiltonian H(t − tn), where H(t) is given in Eq. (2). Specifically, we set the effective Zeeman fields to BJ(t) = B(t − tn)τJ(tn) and the “system-bath couplings” to gJ(t) = g(t − tn), with the functions g(t), B(t) as shown in Fig. 3 and given explicitly in Methods subsection “Smooth evolution of time dependent couplings”.

    Fig. 3: Adiabatic evolution of the time dependent couplings g and B.
    figure 3

    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(tn+1)}. A new cycle then begins from step b.

Importantly, the flux operators Wi,j commute with H(t). Ideally, they therefore retain their initial values, Wi,j = 1. In the presence of decoherence on a real quantum simulator, Wi,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, B0, g1 and T: B0 should be chosen to be large enough compared to the gap Egap 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 g1 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 B0 and g1: 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 gJ = 0, the bath fermions cx and cy are decoupled from the system fermions, cz, and are initialized in the ground state of the bath Zeeman Hamiltonian. Sweeping BJ(t) downwards induces level crossings between states of the system (cz) and bath fermions, and when gJ ≠ 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. 27 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):

$$dot{varepsilon }(t)=-Cleft[varepsilon (t)-{varepsilon }_{{{rm{s}}}}right],$$
(11)

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, B0 and g1, 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. (11), we obtain ε(t) = εs + [ε(0) − εs]eCt. 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 BJ’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, Egap, 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/Egap, 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 B0 = 7, g1 = 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}_{{{bf{J}}}}^{z}) fermions is equal to 1. To probe the convergence to the ground state, we define the quantity νl:

$${nu }^{l}=frac{1}{2pi i}int{{rm{Tr}}}left[{R}^{l}left[{partial }_{{k}_{x}}{R}^{l},{partial }_{{k}_{y}}{R}^{l}right]right]d{k}_{x}d{k}_{y},$$
(12)

where Rl 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 Rbath. 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.

Fig. 4: A single cooling cycle performed on the KSL system.
figure 4

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 g1 = 0.5, B0 = 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 code26 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.

Fig. 5: Cooling in the presence of noise.
figure 5

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−3 errors per cycle per qubit. The cycle duration is set to T = 20 and Nt = 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.

Fig. 6: Steady state in the presence of noise.
figure 6

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/(2n0) taking into account exponentially decaying correlations over approximately n0 = 3 cycles. The cycle duration is set to T = 10 and Nt = 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 wn, and calculate its exponential moving average

$${bar{w}}_{n}=frac{{sum }_{mle n}{w}_{m}{e}^{-frac{n-m}{{n}_{0}}}}{{sum }_{mle n}{e}^{-frac{n-m}{{n}_{0}}}},$$
(13)

where 1/n0 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 n0 = 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 cycle20.

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 uI,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.

Discussion

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. 23).

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, Nt, to optimize the trade-off between diabatic errors and Trotterization noise (which decrease upon increasing T and Nt) versus the noise deriving from environmental decoherence and low fidelity gates (which are typically proportional to the number of gates applied, and hence to Nt). 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 Nt.

Our protocol offers a parametric advantage compared to other recently-proposed methods to prepare chiral 2D states. Ref. 30 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. 31, 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. 23 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. 23 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.

Methods

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 cz. Fluxes can be removed by measuring the plaquette operators Wi,j which correspond to local 6-body spin operators shown in Eq. (10). The excitations Wi,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 code26). 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 decoder32. 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. 27, defined as

$${{mathcal{S}}}(x) = left{begin{array}{ll}0, hfill & x , < , 0 \ 1,hfill & x , > , 1 \ frac{1}{1+{e}^{frac{1}{x}+frac{1}{x-1}}}, hfill &, {rm{otherwise}} ,end{array}right..$$
(14)

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

$$B(t)={B}_{0}left[1-{{mathcal{S}}}left(frac{t}{T}right)right],$$
(15)
$$g(t)={g}_{1}{{mathcal{S}}}left(frac{t}{{t}_{1}}right)left[1-{{mathcal{S}}}left(frac{t-{t}_{2}}{T-{t}_{2}}right)right].$$
(16)

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. 27.

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 uI,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 Nt 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. 27. 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 Jx and to Jy 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 ζ/Nt. 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.

Related Articles

Optical sorting: past, present and future

Optical sorting combines optical tweezers with diverse techniques, including optical spectrum, artificial intelligence (AI) and immunoassay, to endow unprecedented capabilities in particle sorting. In comparison to other methods such as microfluidics, acoustics and electrophoresis, optical sorting offers appreciable advantages in nanoscale precision, high resolution, non-invasiveness, and is becoming increasingly indispensable in fields of biophysics, chemistry, and materials science. This review aims to offer a comprehensive overview of the history, development, and perspectives of various optical sorting techniques, categorised as passive and active sorting methods. To begin, we elucidate the fundamental physics and attributes of both conventional and exotic optical forces. We then explore sorting capabilities of active optical sorting, which fuses optical tweezers with a diversity of techniques, including Raman spectroscopy and machine learning. Afterwards, we reveal the essential roles played by deterministic light fields, configured with lens systems or metasurfaces, in the passive sorting of particles based on their varying sizes and shapes, sorting resolutions and speeds. We conclude with our vision of the most promising and futuristic directions, including AI-facilitated ultrafast and bio-morphology-selective sorting. It can be envisioned that optical sorting will inevitably become a revolutionary tool in scientific research and practical biomedical applications.

Enantioselective synthesis of chiroplasmonic helicoidal nanoparticles by nanoconfinement in chiral dielectric shells

Helicoid metal nanoparticles with intrinsic chirality have unveiled tailorable properties and unlocked many chirality-related applications across various fields. Nevertheless, the existing strategies for enantioselective synthesis of helicoid metal nanoparticles have been predominantly limited to gold. Here, we demonstrate a robust and versatile strategy for the enantioselective synthesis of helicoid nanoparticles beyond gold, leveraging chiral nanoconfinement provided by chiral SiO2 or nanoshells. The chiral nanoconfinement strategy enables the decoupling of ligand-directed crystal growth from chiral induction, allowing for the independent tuning of these two critical aspects. As a result, this approach can not only facilitate the replication of chiral shapes from the chiral nanoshells but also allow the generation of alternative chiral shapes. By employing this approach, we demonstrate the enantioselective synthesis of helicoid Pt, Au@Pt, Au@Pd, Au@Ag, and Au@Cu nanoparticles. The chiroplasmonic properties of Pt- and Pd-based chiral nanoparticles have been discovered, and the inversion of chiroplasmonic properties of Ag-based chiral nanoparticles via facet control has been documented and theoretically explained. The chiral nanoconfinement strategy enriches the toolbox for creating chiral nanoparticles and supports their exploration in diverse applications.

Self-reports map the landscape of task states derived from brain imaging

Psychological states influence our happiness and productivity; however, estimates of their impact have historically been assumed to be limited by the accuracy with which introspection can quantify them. Over the last two decades, studies have shown that introspective descriptions of psychological states correlate with objective indicators of cognition, including task performance and metrics of brain function, using techniques like functional magnetic resonance imaging (fMRI). Such evidence suggests it may be possible to quantify the mapping between self-reports of experience and objective representations of those states (e.g., those inferred from measures of brain activity). Here, we used machine learning to show that self-reported descriptions of experiences across tasks can reliably map the objective landscape of task states derived from brain activity. In our study, 194 participants provided descriptions of their psychological states while performing tasks for which the contribution of different brain systems was available from prior fMRI studies. We used machine learning to combine these reports with descriptions of brain function to form a ‘state-space’ that reliably predicted patterns of brain activity based solely on unseen descriptions of experience (N = 101). Our study demonstrates that introspective reports can share information with the objective task landscape inferred from brain activity.

Dynamic thermalization on noisy quantum hardware

Emulating thermal observables on a digital quantum computer is essential for quantum simulation of many-body physics. However, thermalization typically requires a large system size due to incorporating a thermal bath, whilst limited resources of near-term digital quantum processors allow for simulating relatively small systems. We show that thermal observables and fluctuations may be obtained for a small closed system without a thermal bath. Thermal observables occur upon classically averaging quantum mechanical observables over randomized variants of their time evolution that run independently on a digital quantum processor. Using an IBM quantum computer, we experimentally find thermal occupation probabilities with finite positive and negative temperatures defined by the initial state’s energy. Averaging over random evolutions facilitates error mitigation, with the noise contributing to the temperature in the simulated observables. This result fosters probing the dynamical emergence of equilibrium properties of matter at finite temperatures on noisy intermediate-scale quantum hardware.

Observation of non-Hermitian topological synchronization

Non-Hermitian topology plays a pivotal role in physical science and technology, exerting a profound impact across various scientific disciplines. Recently, the interplay between topological physics and nonlinear synchronization has aroused a great interest, leading to the emergence of an intriguing phenomenon known as topological synchronization, wherein nonlinear oscillators at boundaries synchronize through topological boundary states. To the best of our knowledge, however, this phenomenon has yet to be experimentally validated, and the study of non-Hermitian topological synchronization remains in its infancy. Here, we investigate non-Hermitian topological synchronization, uncovering the influence of system size and boundary site geometry on synchronization effects. We demonstrate that simply varying the lattice size allows transitions between three distinct types of non-Hermitian topological synchronization. Furthermore, we reveal that the geometry of the boundary sites introduces a degree of freedom, enabling the control over the configuration of non-Hermitian topological synchronization. These findings are experimentally validated using non-Hermitian nonlinear topological circuits. This work significantly broadens the scope of nonlinear non-Hermitian topological physics and opens new avenues for the application of synchronization phenomena in future technologies.

Responses

Your email address will not be published. Required fields are marked *