## Abstract

Simulating computationally intractable many-body problems on a quantum simulator holds great potential to deliver insights into physical, chemical, and biological systems. While the implementation of Hamiltonian dynamics within a quantum simulator has already been demonstrated in many experiments, the problem of initialization of quantum simulators to a suitable quantum state has hitherto remained mostly unsolved. Here, we show that already a single dissipatively driven auxiliary particle can efficiently prepare the quantum simulator in a low-energy state of largely arbitrary Hamiltonians. We demonstrate the scalability of our approach and show that it is robust against unwanted sources of decoherence. While our initialization protocol is largely independent of the physical realization of the simulation device, we provide an implementation example for a trapped ion quantum simulator.

## INTRODUCTION

Quantum simulation is an emergent technology that can potentially solve important open problems related to high-temperature superconductivity, interacting quantum field theories, or many-body localization (*1*). While a series of experiments demonstrated the successful implementation of Hamiltonian dynamics within a quantum simulator (*2*–*14*), these works had the simulator initialized in an easily accessible state such as a product state. Consequently, adiabatic evolution from an initial Hamiltonian whose ground state can be prepared to the final Hamiltonian of interest has been used. However, this approach becomes challenging across quantum phase transitions, especially if the transition is of first order.

Our strategy to overcome this problem builds on the recent advances in using dissipative quantum systems to engineer interesting many-body states as the attractor states of such an open quantum many-body system (*15*–*25*). In the past, these dissipative state engineering schemes have been limited to ground states of stabilizer or frustration-free Hamiltonians (*16*, *17*, *26*, *27*), whose ground state can be found by performing local optimizations alone. Unfortunately, almost all many-body Hamiltonians of interest lie outside this class, requiring generalization of the dissipative state preparation procedure.

Here, we present a previously unexplored paradigm for the dissipative initialization of a quantum simulator. We consider a coupling of the many-body system performing the quantum simulation to an auxiliary particle that is dissipatively driven. Crucially, the energy splitting within the auxiliary particle is chosen such that it becomes resonant with the many-body excitation gap of the system of interest, i.e., the difference of the ground-state energy and the energy of the first excited state. Under such a resonance condition, the energy of the quantum simulator is efficiently transferred to the auxiliary particle such that the former is cooled sympathetically (*23*, *25*). Although this setup is only resonant at a single energy, the density of states increases exponentially with energy, resulting in the lowest-lying excitations being the bottleneck for fast ground-state preparation (see the Supplementary Materials for details). While the value of the many-body excitation gap is usually unknown before performing the simulation, we demonstrate that the gap can actually be determined from the quantum simulation data in a spectroscopic measurement. Hence, the dissipative initialization process provides important information about the many-body system of interest at the same time. Notably, we show that the cooling by a single auxiliary particle is efficient, and it is especially robust against unwanted noise processes occurring in the quantum simulator.

To be explicit, we consider different paradigmatic one-dimensional (1D) spin 1/2 many-body systems coupled to a single dissipatively driven auxiliary bath spin (see Fig. 1). This setup can be readily generalized to bosonic or fermionic many-body systems with a larger local Hilbert space, to settings incorporating several bath particles, and to higher spatial dimensions. In the following, we assume a 1D chain of *N* spins governed by the Hamiltonian *H*_{sys}. One boundary spin of the system is coupled to the auxiliary bath spin via an interaction Hamiltonian of the form *g*_{sb} is the strength of the system-bath interaction and the σ* _{i}* refer to the Pauli matrices. The exact values of the dimensionless parameters

*f*are not particularly important. In the models studied here, we find that it is either favorable to choose them roughly equal or have one dominant contribution. In addition, to avoid any symmetries in the interaction preventing the cooling of certain degrees of freedom, it is beneficial to assign slightly different values to them.

_{i}The Hamiltonian of the bath spin *H*_{bath} is given by *H* = *H*_{sys} + *H*_{bath} + *H*_{int} is the total Hamiltonian of the *N* + 1 spin system (*28*).

We would like to stress that such a setup imposes only modest requirements for an experimental implementation, which works equally well for both analog and digital quantum simulators. In particular, we note that our setup does not require control over individual particles of the quantum simulator. In our case, it is sufficient to merely be able to control the bath particle independently of the rest of the system. In addition, the dissipative dynamics can be induced by measuring the spin state of the bath spin followed by a spin flip conditional on measuring the spin in the up state. In Materials and Methods, we give a detailed implementation guide for a trapped ion quantum simulator.

## RESULTS

### Ising chain in a transverse field

As the first paradigmatic model, we consider the Ising model in a transverse field, given by the Hamiltonian*g* is the strength of the transverse field, and *J* is the coupling constant for the Ising interaction. As the Pauli matrices do not commute with each other, it is impossible to minimize the interaction terms and the magnetic field term at the same time, meaning that already this simple model lies outside of the class of frustration-free Hamiltonians. The transverse field Ising model is known to undergo a quantum phase transition at *g* = *J* from a paramagnetic phase (*g* > *J*) to a ferromagnetic phase (*g* < *J*) (*29*). In the following, we will set the energy splitting of the bath spin Δ to be identical to the many-body gap Δ*E* = *E*_{1} − *E*_{0} of the transverse field Ising model, where *E*_{0} (*E*_{1}) is the energy of the ground state (first excited state). In the ferromagnetic phase, the ground state becomes doubly degenerate for large system sizes. Because we are not interested in cooling into a particular ground state, *E*_{1} refers to the first excited state above the ground-state manifold. Below, we will demonstrate that choosing the bath spin splitting as Δ = Δ*E* leads to optimal cooling, and we will show how to extract the (a priori unknown) energy gap Δ*E* from the quantum simulation results.

Let us now analyze the cooling performance of the setup by tracking the system energy 〈*H*_{sys}〉 of the transverse field Ising model in wave-function Monte Carlo simulations of *N* = 5 spins, initially in the experimentally accessible state of all spins pointing up. Figure 2 shows that the energy of the system decreases rapidly and finally approaches a value that is close to the numerically calculated ground-state energy. The cooling performance depends on the choice of the system-bath coupling *g*_{sb} and the dissipation rate γ. In the following, we assume that the time available for the cooling remains fixed. Then, if *g*_{sb} is too small, the cooling dynamics is very slow. On the other hand, if *g*_{sb} is too large, then the system and the bath spin will become strongly entangled, and the cooling performance is reduced. Similarly, if γ is too small, then the cooling is slowed down in the same way, while a too large value of γ will lead to a quantum Zeno suppression of the energy transfer required for the cooling process. Hence, there should be an optimal choice for *g*_{sb} and γ, which leads to a minimum in energy within the available time.

To find this optimal choice, we use a model-independent quantity to measure the cooling performance. For this, we calculate the fidelity of the state of the system with respect to the ground-state manifold of the transverse field Ising model. The fidelity *f* is given by*30*). As the inset of Fig. 2 (A and B) shows, the ground state can be prepared with more than 90% fidelity for the optimal choice of *g*_{sb} = 1.15 *g* and γ = 1.9 *g*.

We can also relate the fidelity *f* to the system energy 〈*H*_{sys}〉. For this, we introduce a dimensionless excitation energy ε, measured in units of the many-body gap Δ*E*, i.e.

In the low-energy limit ε ≪ 1 and assuming that the excitation energy is mostly concentrated in low-energy excitations, ε is related to the fidelity according to ε = 1 − *f*.

We have also checked that our cooling procedure works independently of the choice of *J*/*g*, i.e., both in the ferromagnetic phase and in the paramagnet, as well as independently of the initial state (see the Supplementary Materials for details). Even in the critical regime (*J*/*g* ∼ 1), where the many-body gap is closing, we observe a similar cooling performance. To substantiate this point and also to demonstrate that our cooling protocol is not limited to a particular model, we turn to the especially challenging case of a critical Heisenberg chain in the following section.

### Antiferromagnetic Heisenberg model

As a second paradigmatic quantum many-body model, we investigate the antiferromagnetic Heisenberg chain, given by the system Hamiltonian

This model exhibits an *SU*(*2*) symmetry and serves as the critical point of a Kosterlitz-Thouless transition when the strength of the σ* _{z}*σ

*interaction is varied (*

_{z}*31*). As the many-body gap vanishes in the thermodynamic limit, this model represents a particularly challenging case for our cooling protocol. In addition, the ground state at the critical point is highly entangled (

*32*); hence, we also test the capability of our cooling protocol to prepare entangled quantum many-body states.

The antiferromagnetic Heisenberg model adds one minor complication concerning its ground state preparation compared with the Ising model. Because of an approximate symmetry conserving certain spin-wave excitations, the ground-state cooling performance is limited when the system-bath coupling is restricted to the last spin of the chain. We resolve this issue by an additional system-bath coupling of strength *g*_{sb}/2 to the second last spin of the chain. Figure 3 shows the cooling performance in terms of the system energy 〈*H*_{sys}〉, with an initial state of spins pointing up and down alternately, as a function of the splitting of the bath spin Δ. As in the case of the transverse field Ising model, 〈*H*_{sys}〉 decreases rapidly and reaches a final value that is close to the ground-state energy *E*_{0}. In addition, the cooling is optimal when Δ is chosen to be identical to the many-body gap Δ*E* (*f* = 0.9). Hence, experimentally measuring *H*_{sys} as a function of Δ allows one to obtain the value of the many-body gap Δ*E*, which in itself is an important quantity to understand a quantum many-body system. We also find the final state to be highly entangled (see the Supplementary Materials for details).

However, on many quantum simulation architectures, it might be difficult to experimentally measure the system energy *H*_{sys}, as this will typically require one to perform tomography on all the operators that appear in the system Hamiltonian. Further challenges arise in architectures where not all coupling constants in the Hamiltonian can be perfectly controlled, leading to additional uncertainties in the estimated value of Δ*E*.

Fortunately, it is possible to obtain Δ*E* by measuring only the bath spin. The key idea is to measure the energy *E*_{dis} that is dissipated during the cooling dynamics. Crucially, this energy is related to the number of quantum jumps *N*_{jump} by the relation *E*_{dis} = *N*_{jump}Δ, as a quantum jump will lower the energy of the bath spin by Δ. We note that there are two different ways to obtain *N*_{jump}. First, one can directly count the number of quantum jumps, e.g., by counting the number of emitted photons, if the dissipative flip of the bath spin is realized by a spontaneous emission event. In many setups, however, collecting each emitted photon with high probability might be too challenging. However, as a second method, one can also obtain *N*_{jump} via the integrated probability to find the bath spin in the up state according to*t*_{p} is the total preparation time. As shown in Fig. 3, the minimum of *E*_{dis} is almost identical to the minimum in *H*_{sys}, corresponding to the case where the splitting of the bath spin Δ is identical to the many-body gap Δ*E*. We note that if the system-bath coupling *g*_{sb} or the dissipation rate γ is chosen too large, then the difference between the minima in 〈*H*_{sys}〉 and *E*_{dis} becomes appreciably larger. We also observe that *E*_{dis} is slightly larger in magnitude than the system energy; this can be attributed to the fact that even in the limit of large times, a finite probability for quantum jumps remains, as the ground state of the system Hamiltonian is not a perfect dark state of the quantum master equation (*33*) due to the finite system-bath coupling *g*_{sb}. These additional jumps can also happen for nonoptimal values of Δ, leading to a broadening of the dissipated energy *E*_{dis} in Fig. 3B compared with the system energy 〈*H*_{sys}〉.

### Efficiency of the cooling protocol

For any quantum-state preparation protocol, it is crucial to determine how its properties behave when the size of the system is increased. A protocol is called efficient when the resources required (i.e., the preparation time) grow at most polynomially with the system size. To determine the scaling with system size in an unbiased way, we compute the preparation time *t*_{p} that is required to cool the system down to a fixed dimensionless energy ε, while the system-bath coupling *g*_{sb} and the dissipation rate γ are chosen such that the cooling is optimal. Within our numerical simulations, we use a standard nonlinear optimization scheme (see Materials and Methods for details). In an actual quantum simulator, one can use a hybrid algorithm in which the energies measured on the quantum device are fed back into the classical optimization algorithm (*34*).

Figure 4 shows the scaling behavior of *t*_{p} for the transverse field Ising model. Although the system is cooled across the phase transition into the ferromagnet, the preparation time grows only polynomially with the system size. The same scaling is also observed for the antiferromagnetic Heisenberg model (see the Supplementary Materials for details). This scaling behavior underlines that our cooling procedure is already scalable when using only a single bath spin. As the number of particles is often a scarce resource in a quantum simulator, the required minimal overhead for the initialization allows us to use almost all of the particles for the actual quantum simulation.

### Performance under decoherence

So far, the only source of decoherence in our considerations stems from the dissipative flips of the bath spin. However, in most quantum simulation architectures, there will also be unwanted decoherence processes in the system performing the quantum simulation. Therefore, it is crucial to determine the consequences of this additional decoherence on the performance of our cooling protocol.

As an additional source of decoherence, we consider σ* _{z}* spin flips in the quantum simulation of the transverse field Ising model, applied with a rate κ to all

*N*spins of the quantum simulator. In the ferromagnetic phase, such a spin flip will create two neighboring domain-wall excitations; i.e., when applied to the ground state, the dimensionless energy will approximately increase to ε ≈ 2. This type of decoherence represents a worst case scenario of all local decoherence processes. Hence, we expect that this scenario is quite generic and that our findings should also apply to other many-body models.

To analyze the consequences of these additional decoherence channels, we consider the quantity κ*t*_{p}, which is essentially the probability of any spin to undergo a decoherence event during the preparation time. Then, tracking how the energy ε behaves as a function of κ*t*_{p} allows us to assess the robustness of our cooling protocol under additional decoherence.

Figure 5 shows the system energy for different decoherence rates, from which the behavior of ε is calculated. Crucially, we find that the system contains one excitation, ε ≈ 1 at a value of κ*t*_{p} ≈ 2. This means that the system picks up one excitation when on average all the spins have undergone a decoherence event. This is in stark contrast to the scaling observed in adiabatic state preparation protocols, where the error probability is typically given by the probability that a single spin undergoes a decoherence event, i.e., proportional to *N*κ*t*_{p} (*35*). This improved robustness against decoherence can be attributed to the fact that our state preparation protocol itself is dissipative and therefore can self-correct decoherence events.

### Experimental realization

The proposed initialization protocol can be implemented in a trapped ion system with state-of-the-art technology, e.g., by confining a 1D ion string in a linear Paul trap. Here, we propose an implementation with ^{40}Ca^{+} ions in a setup similar to the one described in (*36*). The spin states are encoded in the optical qubit, ∣↓〉 = ∣*S*_{1/2}, *m* = + 1/2〉 and ∣↑〉 = ∣*D*_{5/2}, *m* = +5/2〉 (see Fig. 1C) with an energy splitting of ℏω_{0}, coherently manipulated by radial laser beams; e.g., the rightmost ion serves as the bath spin (index *b*), while its laser-induced coupling to the neighboring ion (index *s*) implements the system-bath coupling. The bath ion can be isolated from the system interaction by shelving the population to an auxiliary state ∣aux〉* _{b}* = ∣

*D*

_{5/2},

*m*= −5/2〉

*with a laser beam addressing only the bath ion. An experimental realization requires the implementation of the system and system-bath Hamiltonians. For simplicity, we suggest to implement*

_{b}*H*

_{sys}and

*H*

_{sb}in an interleaved fashion by trotterizing the total interaction (

*6*,

*37*).

In trapped ion systems, *H*_{sys} for the transverse field Ising model (*5*) has been realized with up to 53 qubits (*12*). For this purpose, a global bichromatic laser beam with frequency ω_{0} ± δ implements a gate operation by coupling to all radial modes. If δ is larger than the center-of-mass mode frequency, then the resulting spin-spin coupling coefficient shows a power law scaling *J*_{i, j} ∝ 1/∣*i* − *j*∣^{α} (*38*), where α can be varied between 0 and 3 by changing the radial confinement. Implementation of the Heisenberg model is possible by interleaving the spin-spin coupling gates with single-qubit rotations performing a basis change from σ* _{x}* to σ

*and σ*

_{y}*.*

_{z}We propose to implement *H*_{sb} with a separate laser that provides single ion addressing for the bath spin and the neighboring system spin. A Mølmer-Sørensen gate (*39*, *40*) on the radial motional modes bridges two different energy gaps, ω_{s} and ω_{0}, similar to a two-species gate (*41*), and provides a _{0} ± δ, and for the system spin will be ω_{s} ± δ with ω_{s} = Δ*E*/ℏ for optimal cooling. Tuning the latter frequency corresponds to searching for the resonance condition described in the main text. Again, _{sb} ∼ 1. This comes at the expense of an additional interaction between the last and the second-last spins of the system, which is appreciably weaker than *J* and can therefore be neglected. This additional coupling may also be canceled using an additional addressing laser.

Assuming Δ*E* is already known, repumping from ∣↑〉* _{b}* to P

_{3/2}and a subsequent spontaneous decay to ∣↓〉

*on the bath ion can be used to provide a channel for dissipation. The strength of dissipation, γ, within the trotterized scheme can be adjusted by the repumping laser intensity, i.e., the repumping probability during each Trotter cycle. For determination of Δ*

_{b}*E*by recording

*N*

_{jump}, every scattered photon during the repump process has to be detected. This is accomplished by an electron shelving scheme in which the population in ∣↓〉

*is hidden in state ∣aux〉*

_{b}*and a potentially scattered photon bringing the bath ion from ∣↑〉*

_{b}*to ∣↓〉*

_{b}*is detected by measuring fluorescence on the ∣↓〉*

_{b}*(S*

_{b}_{1/2})-to-P

_{1/2}transition. To avoid a perturbation of the system spins, the detection laser has to be tightly focused onto the bath ion.

To be more specific, we assume 5^{40}Ca^{+} ions in a linear chain with single ion axial and radial trapping frequencies of ω* _{z}* = 2π × 0.15 MHz and ω

*= 2π × 0.5 MHz, respectively (*

_{r}*36*). With a resonant Rabi frequency of 2π × 125 kHz for all ions,

*J*

_{i, j}ranges between 2π × 2.7 and 2π × 1.2 kHz for a detuning of δ − ω

_{r}≈ 2π × 10.5 kHz, while the largest Lamb-Dicke parameter is given by η

_{max}= 0.128. For these parameters, the spacing between the bath spin and the nearest system spin of around 14 μm is sufficiently large to provide a factor of 10

^{−7}suppression of the scattering rate for the electron shelving detection on the neighboring ion for a beam focused to 2.6 μm on the bath ion.

In such a setup, the dominant decoherence mechanism is arising from global magnetic field fluctuations. This process can be expressed in terms of a jump operator of the form * _{c}* = 3.3 Hz (

*42*). Figure 6 shows the cooling of such a system to an optimized ground-state fidelity of

*f*= 0.92 in the decoherence-free case, while the presence of decoherence leads to a fidelity of

*f*= 0.89.

An alternative to single ion addressing is to use another isotope for the bath ion, such as ^{44}Ca^{+}. The large isotope shifts of 850 MHz on the S_{1/2}-P_{1/2} transition and 5.3 GHz on the qubit transition (*43*, *44*) will appreciably relax the focusing requirements at the expense of having to achieve an appropriately ordered ion crystal (*45*).

## DISCUSSION

Here, we demonstrated how adding a dissipatively driven auxiliary particle can sympathetically cool a quantum simulator into low-energy states. Our approach is efficient even when using only a single bath spin, and it exhibits strong robustness against unwanted decoherence occurring in the quantum simulator. Future directions include investigating the scaling behavior when optimally varying the coupling constants of the bath in time and when adding multiple bath spins. In the latter case, it will also be of interest to choose different splittings of the bath spins, allowing engineering of tailored bath spectral functions for the quantum simulator.

## MATERIALS AND METHODS

### Numerical simulations

All numerical simulations were performed using a wave-function Monte Carlo approach provided by the QuTiP library (*46*), extended to a massively parallelized version (*47*). Results were obtained by averaging over 1000 Monte Carlo trajectories. We note that we are interested in the long time limit of a weakly dissipative system, i.e., a regime where tensor network algorithms are breaking down (*48*). Numerical optimization of the coupling constants *g*_{sb} and γ was carried out using a Nelder-Mead algorithm. We typically obtain convergence within approximately 50 runs of the simulation, which does not appreciably depend on the size of the system.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/10/eaaw9268/DC1

Section S1. Energy-level representation of the cooling protocol

Section S2. Cooling in the paramagnetic and the critical regime of the Ising model

Section S3. Dependence of the cooling performance on the initial state

Section S4. Efficiency of the cooling protocol for the Heisenberg model

Section S5. Entanglement measure for the ground-state cooling of the Heisenberg model

Fig. S1. Possible paths via which an excitation can be cooled down to the ground state.

Fig. S2. Cooling dynamics in different regimes.

Fig. S3. Cooling performance of the transverse field Ising model in the ferromagnetic phase for various initial states.

Fig. S4. Scalability of the protocol for the antiferromagnetic Heisenberg model.

Fig. S5. Negativity as a measure of entanglement of the prepared states in time for a system of *N* = 6 spins.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was funded by the Volkswagen Foundation and the DFG within SFB 1227 (DQ-mat, projects A01, A04, and B05).

**Author contributions:**M.R. and H.W. designed the cooling protocol. M.R. performed the numerical simulations. F.W., C.O., and P.O.S. designed the experimental implementation proposal. All authors contributed to the writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).