Simulating open quantum systems: from many-body interactions to stabilizer pumping

In a recent experiment, Barreiro et al (2011 Nature 470 486) demonstrated the fundamental building blocks of an open-system quantum simulator with trapped ions. Using up to five ions, dynamics were realized by sequences that combined single- and multi-qubit entangling gate operations with optical pumping. This enabled the implementation of both coherent many-body dynamics and dissipative processes by controlling the coupling of the system to an artificial, suitably tailored environment. This engineering was illustrated by the dissipative preparation of entangled two- and four-qubit states, the simulation of coherent four-body spin interactions and the quantum non-demolition measurement of a multi-qubit stabilizer operator. In this paper, we present the theoretical framework of this gate-based (‘digital’) simulation approach for open-system dynamics with trapped ions. In addition, we discuss how within this simulation approach, minimal instances of spin models of interest in the context of topological quantum computing and condensed matter physics can be realized in state-of-the-art linear ion-trap quantum computing architectures. We outline concrete simulation schemes for Kitaev's toric code Hamiltonian and a recently suggested color code model. The presented simulation protocols can be adapted to scalable and two-dimensional ion-trap architectures, which are currently under development.


Introduction
In view of the inherent difficulties to efficiently simulate quantum physics of an interacting many-body quantum system on a classical computer due to the Hilbert space growing exponentially with the system size, Feynman proposed the idea of a quantum simulator.He suggested a controllable quantum device to efficiently study the dynamics of another quantum system of interest [1].This idea was later refined and formally developed by Lloyd [2] and others, who showed that many-body quantum systems can indeed be simulated efficiently, as long as they evolve according to local interactions.Since then, quantum simulation has become a very active and rapidly evolving research field on its own (see references [3,4] for a recent overview).Driven by remarkable experimental progress and novel theoretical ideas for various physical platforms in recent years, in particular AMO systems ranging from cold atoms [5,6,7,8] and polar molecules [9,10] over trapped ions [11] to photonic setups [12,13,14] and nuclear magnetic resonance [15] have been under investigation for quantum simulation.Similar promising developments have been reported for solid-state systems [16] such as, e.g. , arrays of coupled superconductors [17,18], quantum dots [19] and nitrogen-vacancy centers in diamond [20].
For closed many-body quantum systems, which are well-isolated from their environment, powerful techniques have been developed to control the internal, coherent dynamics.The ability to engineer and tune the underlying single-particle and interaction Hamiltonian terms has enabled the simulation of different classes of quantum many-body models over wide ranges of parameters.Ultimately, though, every quantum system is unevitably coupled also to its surrounding environment.Recently, quantum control of open many-body systems, which amounts to engineering both the Hamiltonian time evolution of the many-body system itself as well as its coupling to the environment [2,21], has become a major research focus.Whereas typically the system-environment coupling leads to detrimental effects on many-body or multi-qubit open systems [22,23,24,25], the ability to control and tailor the associated dissipative processes has been identified as a useful resource: it allows one to dissipatively prepare entangled quantum states and correlated quantum phases from arbitrary initial states [26,27,28,29,30,31], and can also be exploited for dissipative quantum computing [32] and quantum memories [33].
Recently, the elementary building blocks of such an open-system quantum simulator have been shown in an experiment with up to five ions [34].In their work, Barreiro et al. demonstrated the ability to engineer coherent and dissipative multi-qubit quantum operations by the dissipative preparation of Bell states and multi-qubit stabilizer states, the simulation of coherent four-body spin interactions and a quantum non-demolition measurement of four-qubit stabilizer operators.Since the theoretical concepts and details of this work are of general interest to the ion trap community in the context of quantum simulation of spin systems, we provide in the present paper the theoretical framework of the simulation scheme.The present work is motivated by the developments of ideas in the context of topological spin models in the context of quantum computing and condensed matter and the question to what extent these ideas can be realized in existing experimental setups, in particular with linear ion-trap architectures.We focus on the following questions: What are interesting simulation possibilities in state-of-theart ion trap quantum computing setups with moderately large chains of a few, possibly up to a few tens of ions?And how can the currently available experimental resources be exploited in an optimal and experimentally efficient way that allows one to access the physics of minimal instances of complex spin models (as schematically shown in figure 1) with today's technology?Lattice spin models of interest for the gate-based ("digital") quantum simulation with trapped ions.(a) In Kitaev's toric code [35] spins located around vertices of a two-dimensional square lattice interact via four-body interactions ∼ σ x 1 σ x 2 σ x 3 σ x 4 , whereas spins around plaquettes experience z-type interactions, as e.g.∼ σ z 2 σ z 3 σ z 5 σ z 6 .(b) Small instance of a color code spin system, as proposed in [36].Here, spins are located on the sites of a three-colorable lattice interact via four-body plaquette interactions such as ∼ σ (c) Mesoscopic instances of spin models can be mapped onto linear chains of trapped ions, where the spin degree of freedom is encoded in (meta)stable electronic states.Coherent and dissipative time evolution can be simulated by sequences of highly parallel multi-ion Mølmer-Sørensen (MS) gates applied to all (or subsets of) ions, in combination with single-qubit rotations on individual ions and optical pumping of an ancilla ion.
Below, we present a toolbox for the simulation of general Markovian open-system dynamics of mesoscopic spin systems.Our set of tools includes the fundamental building blocks for the simulation of coherent n-body spin interactions, dissipative n-qubit reservoir engineering, and the ability for quantum-non-demolition (QND) measurements of n-particle observables.The simulation scheme strongly makes use of the welldeveloped set of tools for the purpose of quantum state preparation, manipulation and readout of trapped ions [37,11,38].In particular, we show how high-fidelity multiion MS entangling gate operations [39], as first suggested Mølmer and Sørensen, and recently shown for up to 14 ions in the laboratory [40], conveniently bundle the effect of sequences two-qubit operations.This allows one to reduce the experimental simulation complexity significantly and to realize, e.g., coherent n-body interactions in a minimal number of steps.In our simulation architecture, we use optical pumping on individual ions -in combination with coherent gate operations -to tailor the coupling of the spin system to its environment and thereby engineer dissipative n-body quantum processes.
Our "digital" simulation scheme is based on the stroboscopic application of sequences of coherent gate operations in combination with dissipative time steps to realize open-system dynamics.It complements existing proposals of quantum simulation with ground state ions [41,42,43] or ions excited to Rydberg states [44].In these "analog" quantum simulators, the common principle is to use externally controllable fields to engineer effective "always-on" Hamiltonians, which microscopically realize the model of interest directly.Recently, remarkable first experiments have demonstrated the simulation of (relativistic) single-particle dynamics in an external potential [45,46,47] and experimental studies of the physics of few interacting Ising spins [48] under frustration [49].
We point out that the presented "digital" simulation scheme is suited for the simulation of mesoscopic spin systems corresponding to chains of up to a few tens of ions, which with state-of-the art ion trap technology can be controlled accurately.However, similar protocols can be realized in scalable and two-dimensional ion-trap architectures, to whose development currently a lot of effort is devoted [50,38,51,52,53], and also on other physical simulation platforms.In fact, in previous work a "digital" quantum simulation architecture for open-system dynamics of many-body spin models has been developed for neutral Rydberg atoms in optical lattices [54].
In section 2 we introduce the general idea of our simulation architecture and give a concise summary of the main results.The details of the simulation of coherent and dissipative many-body interactions are provided in sections 3 and 4. In section 6 we briefly discuss the effect of imperfections in the simulation scheme.We illustrate our simulation scheme in section 5 for two examples of interest in the context of topological quantum computing, namely small-scale implementations of Kitaev's toric code [35] and a minimal instance of a color-code model [36].We conclude with an outlook.

Open-system dynamics
In the following we are interested in open-system dynamics of many-body quantum systems.The dynamics of an open quantum system which is coupled to an environment can be described by a completely positive Kraus map [55] ρ where ρ denotes the reduced density operator of the system, {E k } is a set of operation elements satisfying k E † k E k = 1, and we assume an initially uncorrelated system and environment.For the case of a closed system, decoupled from the environment, the map (1) reduces to ρ → U ρU † with U the unitary time evolution operator of the system.
In the literature on quantum control of open quantum systems, the required set of operations to realize different classes of quantum operations (1) as well as efficiency and universality aspects have been discussed [21,56].In reference [34], several specific examples of Kraus maps, whose dissipative dynamics can be used for dissipative quantum state preparation, e.g. for pumping into entangled states, have been discussed and implemented experimentally.
The Markovian limit of the general quantum operation (1) for the coherent and dissipative dynamics of a many-particle spin system is given by a many-body master equation for the density operator ρ(t) of the many-body system.The coherent part of the dynamics is described by It is generated from a Hamiltonian H = α H α which is a sum of terms H α , which can in general involve higher order n-body interactions, which act on a quasi-local subset of particles ‡.Dissipative time evolution is described by a Liouvillian part of the master equation The individual terms L diss α ρ are of Lindblad form [57], and are determined by quantum jump operators c α , which either act on single or on subsets of particles, and by respective rates γ α at which these jump processes occur.Engineering open-system dynamics thus amounts to designing and engineering couplings of the quantum system to its environment, such that the resulting many-particle dynamics is then governed by discrete Kraus maps or master equations with quasi-local Hamiltonian and dissipative terms.

Many-body quantum systems: Kitaev's toric code as a representative example
In the following, we will in consider the simulation of many-body lattice spin models, which are of interest in the context of topological quantum computing and memories.As a paradigmatic example of this class of spin models we discuss in some detail Kitaev's toric code Hamiltonian, which is sketched in figure 1(a).This model exemplifies in a transparent way the challenges that one encounters also in the quantum simulation of ‡ Throughout this article we set = 1.
related models, such as e.g. in a recently suggested color code model (see figure 1(b)), which we discuss in more detail in section 5.2.
In Kitaev's toric code model, as sketched in figure 1(a), spins are located on the edges of a two-dimensional square lattice.The Hamiltonian is given by H = −E( s A s + p B p ), which is a sum of stabilizer operators which describe four-body interactions of spins, which are located around the vertices (stars) s and plaquettes p of the lattice, respectively.All four-body stabilizers have eigenvalues ±1 and mutually commute.The ground state(s) is/are thus given by the simultaneous eigenstate(s) of all stabilizers with eigenvalues +1 (assuming E s , E p > 0).The degeneracy of the ground state depends on the boundary conditions and topology of the setup.Excited states in this model correspond to violations of these stabilizer constraints, i.e., −1 eigenstates with respect to either the A s or B p stabilizers.They can be associated with localized quasiparticles residing on the corresponding vertices and plaquettes of the lattice (as illustrated in figure 5(b)).They exhibit anyonic statistics under braiding, i.e. when trajectories of different types of quasiparticles are winded around one another.Preparation of the system in the ground state manifold, starting from an arbitrary initial (excited) state, can be achieved by a dissipative dynamics which is governed by a many-body master equation (2) with quantum jump operators  These collective operators act on four spins located around a vertex (site) of the lattice, as depicted in figure 1(a).The index i denotes one arbitrary spin of the four involved spins.A four-body jump operator c α induces dissipative dynamics, which pumps the four spins from the +1 into the -1 eigenspace of A α = σ x 1 σ x 2 σ x 3 σ x 4 (see figure 2).The projector ) applied to any +1 eigenstate of A α vanishes; as a consequence all +1 eigenstates are "dark states" and remain unaffected.In contrast, the spin flip σ z i applied to one of the four spins (e.g.i = 1) can incoherently convert -1 into +1 eigenstates, e.g., Here, |± are the eigenstates of σ x : σ x |± = ±|± .
The above example illustrates that the difficulty to be overcome in simulating the coherent Hamiltonian dynamics lies in finding a way to realize the four-body Hamiltonian interaction terms.The realization of the dissipative "cooling" dynamics into the ground state(s) by means of the described collective dissipative processes requires the engineering of a coupling of the spin system to an artificial, tailored environment.An analog simulation of these coherent and dissipative higher-order n-body interactions, i.e., by a direct engineering using "always-on" external fields, is demanding because these higher-order effective interactions must be constructed from underlying one-and two-body interactions.In this scenarion, typically, the interaction strengths and dissipative rates of the n-body processes, which typically arise in a perturbative limit, are much smaller than dominant one-and two-body interactions.
Therefore, we aim to realize the coherent and dissipative dynamics according to (1) or ( 2) in a digital simulation, i.e. by stroboscopic sequences of gates and dissipative operations.Here, higher-order n-body interactions can be obtained non-perturbatively as leading-order terms from the application of one-, two-or n-body quantum gates.The corresponding interaction strengths are virtually independent of the order n of the interaction terms and ultimately only limited by the gate durations in the underlying quantum circuits.
In the case of continuous time dynamics, we apply these operations over small time steps τ , such that the master equation ( 2) emerges as an effective, coarse-grained description of the time evolution.For small time steps, the time evolution can be implemented through a Trotter expansion of the propagator corresponding to Eq. ( 2) Errors from possible non-commutativity of the quasi-local terms in L are bounded [55] and can be reduced by resorting to shorter time steps τ or higher-order Trotter expansions [58].On the other hand, as we will discuss below, it is also possible to engineer sequences of discrete Kraus maps (1), which can for instance be employed for dissipative quantum state preparation in a minimal number of steps.

Experimental tools for digital quantum simulation with trapped ions
Motivated by the present availability of well-developed set of coherent and dissipative tools [34] in state-of-the-art linear ion-trap architectures [37], we consider a setup in which the spins of a (possibly two-or three-dimensional) lattice model with a mesoscopic number of particles are mapped onto a linear chain of ions, where the spin degrees of freedom are encoded in two (meta-)stable internal states of the ions.Although our approach can be realized with any universal set of gate operations, we focus on a realization, which benefits from highly parallel multi-ion MS gates as the principal building block for the implementation of unitary and dissipative simulation time steps in eq. ( 7).The MS gate operation [39] is based on pairwise two-ion interaction terms (as illustrated in figure 3), and can be parametrized by two angles θ and φ, The sum in the collective spin operators S x,y = n i=0 σ x,y i with σ x,y i the usual Pauli matrices, is understood to be performed over all ions involved in the gate.This multi-ion entangling gate operation is complemented by (non-entangling) single-and multi-qubit rotations, whose physical implementation is discussed, e.g., in [34].In addition to this universal set of coherent gate operations, the use of optical pumping on individual ions (as demonstrated e.g.[59]) constitutes the dissipative ingredient for the engineering of dissipative many-body spin dynamics.(8).All pairs {i, j} of ions involved in the gate interact with equal strength (represented as links).(b) A (4 + 1) ion entangling MS gate applied to 4 system ions and an ancilla ion (index zero) can be used to coherently map the information about whether the 4 system ions are in a +1 (-1) eigenstate of the 4-body interaction term ∼ σ x 1 σ x 2 σ x 3 σ x 4 , onto the logical states |0 and |1 of the ancilla ion.In the Bloch sphere representation this mapping can be understood as a rotation of the ancilla qubit initially prepared in |0 around the x-axis.The rotation angle depends on the state of the system ions, and is chosen such that for any +1 eigenstate, as e.g.| + + + + , the ancilla qubit ends in |0 after the MS gate, whereas for -1 eigenstates such as e.g.| + + + − it is transferred to |1 .This mapping mechanism not only works for 4-body interactions, but can be used for general n.For n odd, the ancilla qubit is transferred to σ y eigenstates.
We have summarized the the basic idea of the simulation of coherent and dissipative dynamics corresponding to Kitaev's code model in figure 4, to be explained in more detail in the following sections.
We will show in more detail in section 3 that the unitary propagators e L coh α τ ρ corresponding to n-body interaction Hamiltonians H α (such as, e.g., the four-body term in Eq. ( 5)) can be implemented efficiently in an experiment (i.e. by a minimal number of gates) by combining standard single-qubit gates with (n + 1)-ion MS gates, which are applied to the n system ions and an additional ion, which encodes an ancilla qubit.Dissipative dynamics according to propagators e L diss α τ ρ with many-body jump-operators c α can be achieved by combining the coherent gate operations with a dissipative step in the form of optical pumping of the ancilla ion.
for a time step τ is efficiently achieved in three steps: (i) First, an entangling MS gate U MS (π/2, 0) applied to the four system ions and the ancilla ion coherently maps the information on whether the system ions are in a +1 or -1 eigenstate of H α onto the logical states |0 and |1 of the ancilla qubit (cf.figure 3).(ii) Second, a single-qubit gate exp(iφσ z 0 ) applied to the ancilla ion effectively imprints a phase −φ (φ) on all +1 (-1) eigenstates of H α .(iii) Finally, the initial mapping is reversed by an inverse MS gate, which disentangles the ancilla from the system ions, which have evolved according to exp(iφσ x 1 σ x 2 σ x 3 σ x 4 ).(b) Dissipative evolution, i.e. "cooling" of the system ions into the +1 eigenspace of H α mediated by Lindblad dynamics with four-body quantum jump operators i ] is applied between the mappings (i) and (iii).In step (ii), all +1 eigenstates of H α are left invariant, whereas a spin flip σ y i with a θ-dependent amplitude is applied to one of the four system ions can convert -1 into +1 eigenstates.The angle θ controls the conversion probability and allows one to tune from small probabilities (θ 1, master equation limit) to unit pumping probability (θ = π/2).After the three steps, generally, the ancilla qubit is entangled with the system ions.(iv) Finally, the ancilla qubit is incoherently reset to its initial state |0 by optical pumping.This dissipative step enables to carry away entropy and "cool" the system qubits.

Simulation of coherent n-body interactions
In this section we describe in more detail the stroboscopic simulation of coherent dynamics e L coh α τ according to n-body spin interaction Hamiltonians.

The Mølmer-Sørensen multi-ion entangling gate
The main resource in our simulation scheme are multi-ion MS gate operations (8), which rely on the application of a bichromatic laser field to the ions [39].The two frequency components are chosen close to the qubit transition, fine-tuned such that an effective second-order coupling between pairs of ions is generated by off-resonantly coupling to the blue and red motional side-bands of the common vibrational center-of-mass mode of the ion string.Within the Lamb-Dicke regime, where the ions are spatially confined to a region much smaller than the wavelength of the qubit transition, the MS gate operation is particularly robust and works in principle without the necessity of cooling to the motional ground state [60].The gate has been successfully demonstrated [61] with remarkably high fidelities (99.3 % for a pair of ions [62]) and recently for strings of up to 14 ions [40].A detailed discussion of the MS gate, in particular regarding the experimental implementation and optimization, can be found in [60].The properties of the MS gate (8) applied to n + 1 ions, which we will repeatedly use in the following, are: • The phase θ is the main control parameter of the gate; for θ = π/2, the gate is maximally entangling, i.e., the computational basis states are mapped to states, which are up to local rotations equivalent to GHZ states [39].Shifting the optical phase of the bichromatic driving field allows one to switch between a σ x -type (φ = 0) and a σ y -type (φ = π/2) MS gate.
• "Backward" MS gates (i.e. for negative values of θ) can be realized by "forward" gates, since for n even, with σj = cos φσ x j + sin φσ y j .In particular, fully entangling MS gates are (up to local rotations for an odd number of ions) equivalent to their inverse operations.Using only "forward" MS gates can be experimentally convenient as a sign change of θ generally requires frequency changes of the driving field.

Circuit decomposition for four-body spin interactions
Let us now outline the procedure to simulate a coherent time step e L coh α τ for a fourspin interaction term H α = −E α A α with A α = σ x 1 σ x 2 σ x 3 σ x 4 (see figure 1).Although the unitary propagator can in principle be implemented with a standard universal set of single-and two-qubit gates available for ions [37], here we use an alternative technique, which harvests the multi-ion MS gates and makes use of an ancilla qubit [63] encoded in an additional ion (see figure 4).This technique has been used in [34] to experimentally realize four-body spin interactions.
The approach consists of a sequence of three gate operations: (i) First, a fully entangling MS gate U MS (π/2, 0), applied to the four system ions and the ancilla ion, coherently maps the information, whether the system ions are in a +1 or -1 eigenstate of A x onto the ancilla qubit (see figure 4(a)).(ii) Second, a single-qubit gate U anc (φ) = exp(iφσ z 0 ) is carried out on the ancilla ion.Due to the previous mapping, this operation on the ancilla qubit is equivalent to manipulations on the +1 and -1 subspaces of A α .(iii) Finally, the mapping is reversed by an inverse MS gate U MS (−π/2, 0) on all ions.The evolution according to the three unitaries is given by with the operator Ŝx = n i=1 σ x i acting on the n system ions.Using the identities and one finds that for n = 4 eq.( 10) indeed reduces to As a consequence, the ancilla -initially prepared in |0 -factorizes out from the dynamics of the system ions, which in turn evolve according to the unitary time evolution operator exp(iφA α ).Here, from exp(iφA α ) = exp(−i(−E α A α )τ ) one identifies the energy scale of the four-body interaction as E α = φ/τ , where τ is the physical time, which is needed to perform all gates required for one full simulation time step (7).Note that pairwise interactions among the system ions, present in the two-body Hamiltonian underlying the MS gate (8), cancel out in the inverse mapping step (second MS gate).

Toolbox for simulation of n-body spin interactions
The simulation scheme outlined above for four-body interactions is readily generalized to arbitrary n-body interactions of the form A = n i=1 σ α i with σ α i ∈ {1, σ x i , σ y i , σ z i }.This is possible by applying local rotations to (a subset of) system ions before and after the gate sequence and thereby effectively transforming σ x i into σ y i or σ z i , and by varying the phase φ of the MS gates.This enables, e.g., the simulation of z-type four-body plaquette interaction terms as required for the toric code Hamiltonian (see figure 1 and section 5.1).The required gate sequences are summarized in table 1.If certain ions are supposed not to participate in the interactions (i.e.σ α i = 1) this can be achieved in different ways: (i) by focusing the driving laser of the MS gate operation only onto the relevant subset of ions, or (ii) by hiding the electronic population of these ions in uncoupled electronic states for the duration of the gate sequence [64], or (iii) by means of refocusing techniques [65]: As shown in reference [66] interspersing MS gates applied to all ions with single-qubit gates on individual ions allows one to decouple effectively certain ions from the dynamics.A set of convenient gate sequences for this purpose is discussed in Appendix B. Circuit decompositions for the simulation of more complex many-spin interactions going beyond n-qubit Pauli operators, such as e.g.ring-exchange interactions, can be worked out and implemented in an analogous way (e.g. in [54] such cases are discussed).
We note that in the gate-based "digital" simulation scheme the energy scale E 0 of the n-body interactions is essentially independent of the order n, and mainly limited by the inverse time required for performing the (n + 1)-ion MS gates.This is in contrast to analog simulation approaches, where higher-order interactions typically arise in a perturbative limit from a two-body Hamiltonian, thus with correspondingly smaller energy scales.
Table 1. Circuit decompositions for the simulation of one time step of coherent dynamics according to the time evolution operator U = exp(iφA).The unitary block is implemented by two MS gates applied to the n system ions and an ancilla qubit (# 0) initially prepared in |0 , and a single-qubit rotation on the ancilla qubit.
Finally, we remark that the coherent n-body interactions A α = σ x 1 . . .σ x n can also be achieved without an ancilla qubit by a slight modification of the employed quantum circuit (see appendix Appendix A for details).

Engineering dissipative many-body dynamics
In this section we show how to engineer dissipative dynamics according to n-qubit stabilizer pumping.To be specific, we first discuss the implementation of master equation dynamics governed by four-body plaquette quantum jump operators, c α = ), as required e.g. for the ground state cooling of Kitaev's toric code (as discussed above in section 2).The stabilizer pumping, as described in this section, has been demonstrated in an experiment with five ions [34], four of them representing four system spins, which can be regarded as one plaquette, and one additional ion encoding an ancilla qubit, which has been optically pumped to engineer the dissipative four-spin dynamics.

Engineering four-body quantum jump operators for stabilizer pumping
The dissipative pumping dynamics to "cool" into the ground state manifold of Kitaev's toric code Hamiltonian, as sketched in figure 2, is implemented by three unitary gate operations applied to the four system ions and the ancilla qubit initially prepared in |0 , followed by optical pumping of the ancilla qubit.The sequence of unitaries is with the two-qubit gate (i) As for the coherent simulation, an entangling MS gate U MS (π/2, 0) first maps the information on whether the four system ions are in the +1 or -1 eigenspace of A α onto the logical states of the ancilla qubit.(ii) Next, the gate C(θ) realizes a spin flip with a θ-dependent amplitude, provided the ancilla is in |1 , i.e.only if the system spins are in a -1 eigenstate of A α .In Appendix C we give a possible decomposition of C i (θ) into global MS gates and single-ion rotations.(iii) After reversing the initial mapping (i) by another (inverse) MS gate, the ancilla qubit is in general entangled with the four system spins.(iv) Finally, optical pumping of the ancilla ion back to its initial state |0 constitutes the dissipative element in the sequence, which renders the dynamics of the four system spins irreversible and enables to carry away entropy and thereby "cool" the system qubits.The unitary sequence ( 14) can be expressed as Here we have made use of the fact that all pairwise interaction terms not involving either the ancilla ion or the i-th system ion cancel out due the inverse MS gate.The resulting operation , where U 2 = cos((π/2) Ŝx )σ z 0 + sin((π/2) Ŝx )σ y 0 and U 3 = U † 0 exp(iθσ y i )U 0 , can be further simplified using the operator identities ( 11) and (12) for n = 4, as well as In combination with the subsequent optical pumping of the ancilla ion, the resulting dynamics is given by the quantum operation with the operation elements With a probability p = sin 2 θ states in the -1 eigenspace of A α are dissipatively converted into +1 eigenstates, while the +1 eigenspace is left invariant by the operation.For θ = π/2 cooling occurs with unit probability.For small values θ π/2 equation ( 17) can be expanded up to second order in θ, yielding the standard form of a Lindblad master equation ( 4) with a four-body jump operator c α and the corresponding dissipative rate γ α = θ 2 /τ .Here, as above, τ is the physical time needed for the implementation of one simulation time step (7).

Toolbox for dissipative quantum simulation
The described four-step scheme is readily generalized to n-body stabilizer cooling with n-qubit quantum jump operators of the form c = 1 2 σ z i (1 − A α ), where A α = n j=1 σ α i with σ α i ∈ {1, σ x i , σ y i , σ z i }.In table 2 the required gate operations and the resulting n-body quantum jump operators are listed.By combining the outlined scheme with local rotations on (subsets of) the system ions, this allows one, e.g., to engineer cooling dynamics according to z-type four-body quantum jump operators ), which are required for ground state preparation in Kitaev's toric code model, as explained in section 5.1.
Required gate operations for the simulation of dissipative dynamics according to n-body quantum jump operators c, which generate pumping into the +1 eigenspace of the many-body stabilizer operator A = σ x 1 . . .σ x n (and σ y 1 . . .σ y n , respectively).The form of the two-qubit gates C(θ) are listed for different values of n; see Appendix C for convenient decompositions into single-qubit and collective MS gate operations.Here |y ± denote the eigenstates of σ y , σ y |y ± = ±|y ± .

Applications
In this section we discuss two examples of possible realizations of complex spin models within the presented simulation scheme.We start with a few comments on boundary effects related to the simulation of minimal instances of Kitaev's toric code Hamiltonian [35].Subsequently, we proceed with the discussion of a minimal instance of a recently suggested color code model [36], whose implementation is feasible in a setup of eight ions.The simulation of other models, typically involving similar many-body spin interaction terms, poses comparable demands with regard to the simulation abilities of coherent and dissipative dynamics.

Kitaev's toric code
For small instances of the toric code model, as e.g. a system of two plaquettes as illustrated in figure 5(a), it is possible to define reduced two-or three-body stabilizers for the spins located at the border of the system.For this small system of two plaquettes, coherent and dissipative dynamics can be implemented in a setup of of 7+1 ions.Dissipative dynamics that cools into the ground state manifold of the model, is realized by a Liouvillian with four-body quantum-jump operators as given in eq. ( 5).As discussed above, they realize pumping into the low-energy subspaces, i.e. eigenspaces corresponding to an eigenvalue +1 of the plaquette and vertex terms of the Hamiltonian.In figure 5(c) the effect of a quantum jump induced by the operator ) is illustrated.If an excitation on the left plaquette is present, in a quantum jump the state of the system spins is of the converted into state, which is a +1 eigenstate of the left plaquette operator.In this process the excitation hops over to the neighboring plaquette.If a second excitation was present on the second plaquette, the pair of excitations is "annihilated" in such step.This removal of a pair of excitation constitutes a cooling event.Note that the other type of excitation remains unaffected by this dynamics.Alternatively, cooling takes also place if one of the spins which forms the border of the system is flipped (e.g.spin # 2 in the lower part of figure 5(c)), since then a single excitation is "pushed out" from the lattice system.For the corners and edges of a small lattice system two-and three-body jump operators for x-type stabilizer pumping can be realized in analogy.

Simulation of a color code model
The idea of storing and processing quantum information in naturally protected quantum systems has attracted a lot of interest in recent years [35,67].Here, protection from local errors is achieved by encoding quantum information not in individual physical qubits, but instead in ground states of topologically ordered quantum systems, which provide an energy gap to excited states and exhibit a ground state degeneracy, which cannot be lifted local perturbations.One class of topological quantum error correcting codes that have been proposed in this context, are color codes [36], which exhibit remarkable computational and error correcting capabilities.In particular, they allow us to implement the Clifford group in a fully topological way within the ground state manifold, without the need of individual addressing of physical qubits or braiding of quasiparticles.
Here, we outline how a minimal instance of such a color code can be realized with the discussed simulation tools, and discuss the state preparation and implementation gate operations, as well as readout.For a detailed introduction to color code models we refer the reader to [68].
A minimal version of a color code, consisting of seven physical qubit located at the corners of three plaquettes, is shown in figure 1(b).Including one ancilla qubit for the implementation of coherent and dissipative dynamics this system could be simulated with a string of eight ions.Qubits located around plaquettes interact via four-body x and z-interaction terms: The Hamiltonian , and similar interaction terms for the other two plaquettes, consists of six mutually commuting stabilizer operators.
Coherent dynamics according to this Hamiltonian can be realized by implementing unitary time steps as outlined in section 3 with the help of an ancilla qubit.Cooling into the ground state manifold can be achieved by a Liouvillian dynamics associated with a set of six four-body jump operators, such as c = , thereby driving the system spins dissipatively into the +1 eigenspaces of the six stabilizers.Alternatively, ground state cooling can be realized by a sequence of deterministic cooling steps: Starting with the system spins in the fully polarized state |0 ⊗7 (being already a +1 eigenstate of all z-stabilizers) it suffices to implement three dissipative Kraus maps (1) such as with A = σ x 1 σ x 2 σ x 3 σ x 4 and accordingly for the other two plaquettes, in order to prepare the system also in the +1 subspaces of the x-type stabilizer operators.As discussed in section 4 this is achieved by choosing two-qubit correcting gates C i (θ = π/2) (see equation ( 15)).
Excited states |ψ of the Hamiltonian H cc correspond to states, where the system spins are in -1 eigenstate(s) with respect to certain stabilizer(s).The quasiparticles of x (or z) type associated with these violations of the stabilizer constraints are located on the corresponding plaquettes, for instance on the uppermost plaquette, if Since there are only six stabilizer constraints for seven system spins, the ground state is degenerate and thus offers the possibility to encode one logical qubit.An operator basis for this logical qubit can be constructed by the global operators X = 7 i=1 σ x i and Ẑ = 7 i=1 σ z i .These two logical operators commute with all six stabilizers of the code, thus they leave the system within the ground state manifold.
The logical qubit can be initialized in the logical state | 0 by (dissipatively) preparing the system -in analogy with the four-body stabilizer cooling -in a +1 eigenstate of the global operator Ẑ, such that Ẑ| 0 = | 0 .The logical | 1 -state is then obtained by the application of the logical X-operator, | 1 = X| 0 , which corresponds to a single-qubit rotation applied to all seven system ions.This minimal color code setup also allows to implement single-qubit gates belonging to the Clifford group in a topological way: The Hadamard H and phase-shift gate K can be implemented by applying the corresponding operations globally to all seven system ions, i.e.Ĥ = 7 i=1 H i and K = 7 i=1 K i .The logical operators then directly fulfill the required transformation properties, as for example Ĥ † X Ĥ = Ẑ.
Remarkably, once the system is prepared in the code space, the logical singlequbit gates can be performed by simple global single-qubit rotations without the need of addressing individual ions.These operators take ground states to ground states, the system stays in the code space, and thus the quantum gates are achieved without braiding of quasiparticles.Similarly, readout measurements can be performed globally, i.e. by standard fluorescence imaging of all ions measured either in the x or in the z-basis.
For the realization of a topological C-NOT gate operation, the minimal system consists of two seven-qubit-layers encoding two logical qubits.Its implementation therefore requires an experimental setup of fifteen ions, and might thus become experimentally feasible in the near future.

Noise and Imperfections
In a stroboscopic simulation of a many-body master-equation ( 2) several sources of imperfections occur.First of all, Trotter errors from the non-commutativity of coherent (and also dissipative) terms arise for each time step of the simulation.These are bounded and can be reduced by resorting to smaller time steps and higher-order Trotter decompositions.In addition, imperfect gate operations in the quantum circuits lead to errors.Their effects on gate-based quantum simulations has discussed in detail e.g. in reference [63].

Generic effect of gate imperfections on the quantum simulation
Here, we first briefly discuss the generic effect by considering a particularly transparent example of a pulse length error in the simulation of coherent dynamics U = exp(iφA) according to a four-body spin interaction A = σ x 1 σ x 2 σ x 3 σ x 4 , as explained in section 3. We assume that the only error is a pulse length error in the single-qubit gate U anc (φ) = exp(iφσ z 0 ) applied to the ancilla qubit in the three-step sequence eq.( 10).In one small time step (φ π/2), the system spins evolve according to Assuming that the actual value φ fluctuates (e.g.due to laser intensity fluctuations) in the experiment from time step to time step according to a Gaussian distribution p(φ) = 1/ √ 2πσ 2 exp[−(φ−φ 0 ) 2 /(2σ 2 )] around the mean value φ 0 with a variance σ φ 0 we obtain, after averaging over φ, the modified equation of motion Thus, one finds dynamics according to a four-body Hamiltonian H eff = −(φ 0 /τ )A, where a systematic shift in φ results in a systematically larger or smaller energy scale.In addition, the stochastic Gaussian fluctuations in φ cause a collective dephasing dynamics (in the σ x -basis), described by a Liouvillian with a dephasing rate γ = σ 2 /τ and a hermitian four-body quantum jump operator A (see last term in (23)).The effect of other gate errors in the circuit decompositions for coherent and dissipative dynamics can be analyzed in an analogous way.A more specific error analysis, going beyond these quite general arguments, requires more precise information about the dominant error sources in concrete experimental setups.

Comparison with experimental stabilizer pumping
In the work [34] four-qubit stabilizer pumping and the effect of errors have been studied experimentally.For the benefit of the reader and to make the present discussion selfcontained we find it worthwhile to review briefly the main findings, as explained in detail in the supplementary information of the Nature article [34], to relate this to the present discussion.In the experiment with five ions (which encoded four system qubits and one additional ancilla qubit) stabilizer pumping with 100% pumping probability per step, from the -1 into the +1 eigenspace of the four-qubit stabilizer operator A = σ x 1 σ x 2 σ x 3 σ x 4 has been applied repetitively.The corresponding discrete Kraus map reads ρ → E 1 ρE † 1 + E 2 ρE † 2 with operation elements Starting with the four system qubits in the initial state |1111 , ideally these reach the four-qubit GHZ state (|0000 + |1111 )/ √ 2 after a single application of the above Kraus map.This is reflected by the fact that the expectation value of the four-qubit stabilizer A assumes a value of +1 after the application of this dissipative step.At the same time, the expectation values of the two-qubit stabilizer operators σ z i σ z j , as depicted in schematically in figure 6, should ideally remain unaffected by the four-qubit stabilizer pumping dynamics and stay at a value +1.
x Figure 6.a) Results from a numerical simulation, which accounts for the effect of single-qubit gate errors in repeated stabilizer pumping according to the dissipative map (24).The stabilizer expectation values are obtained by averaging over 10000 realizations, assuming un-correlated errors from gate to gate, in the phases of the single-qubit rotations exp(−i(θ + δθ)/2σ z 0,4 ) with δθ obeying a Gaussian distribution with zero mean and a variance of δθ of 0.3 × π/2.b) Experimental repeated stabilizer pumping.Plot reproduced from data given in the supplementary information of the Nature publication [34].Quantitative differences from the numerical simulation are mainly due to additional errors in the global gate operations, whose precise form is unknown and which have not been taken into account in the theoretical error model.
In the experiment, the Kraus map (24) has been realized by a quantum circuit consisting of global rotations and MS gate operations applied to all five ions and addressed z-type single-qubit rotations, which only involve the ancilla qubit (index # 0) and the system qubit with index # 4 (see the supplementary information of [34] for the exact form of the experimental circuit decomposition).It is reasonable to assume that errors in the global gate operations affect all ions to a similar degree, whereas gate errors in the addressed single-qubit gates lead to additional errors, which dominantly act on the system qubit # 4. As a consequence, one expects that under a repeated application of the dissipative map (24), the expectation values of the two-qubit stabilizer operators, which involve σ z 4 to decay faster than those not involving the ion #4.This decay behavior, as found from a numerical simulation (see figure 6(a) and figure caption for details), has qualitatively been observed in the experiment [34].In addition to the errors in the single-qubit gates (which experimentally are rooted in pulse length errors and / or laser intensity fluctuations, as discussed in the previous subsection), the dominant source of errors are imperfections in the MS gate operations, which result in quantitative differences in the observed decay of stabilizer expecation values.In the language of stabilizer models, these errors in the simulation correspond to unwanted heating processes with respect to the z-type stabilizers during the four-body stabilizer pumping according to the map (24).

Conclusions and Outlook
In this work we have discussed a toolbox for "digital" quantum simulation with linear chains of trapped ions.We have outlined the theoretical concepts and details of the experiment, which recently demonstrated the building blocks of an open-system quantum simulator with up to five ions.Furthermore, we have discussed how our scheme allows one to explore the physics and simulate the coherent and dissipative dynamics of minimal instances of spin models involving n-body interactions and constraints, such as e.g.Kitaev's toric code and a minimal version of topological color code model.Similarly, circuit implementations for more complex coherent and dissipative n-body interaction terms as, e.g., plaquette exchange interactions can be developed; see for instance [54].
Here, we have focused on open-loop dynamics, where coherent and dissipative time evolution in stabilizer models is implemented with the aid of an ancilla qubit, which is not observed.It is known that such open-loop dynamics involving a single, nonobserved ancilla qubit is not sufficient to realize the most general Markovian multiqubit open-system dynamics.As shown in [21], this can be achieved by a closed-loop simulation scenario.Here, general open-system quantum operations are realized by consecutive sequences of coherent operations applied to the system qubits and the ancilla, interspersed with measurements of the ancilla qubit in an appropriate basis.The gathered information from the outcomes of the sequential ancilla measurements can be classically processed and used for feedback operations on the system.We note that the described scheme also allows one to extract such information about the system qubits via a measurement of the ancilla qubit, as schematically shown in figure 7 via a measurement of the ancilla qubit.The circuit for coherent simulation of a four-body spin interactions as discussed in section 3 (cf.eq. ( 13)) realizes the unitary operation U = exp(iπ/4σ z 0 ⊗ A).Thereby, the ancilla qubit initially prepared in |+ is coherently mapped onto the two σ y -eigenstates, depending on whether the four system qubits are in a +1 or -1 eigenstate of A. This information obtained from a subsequent measurement of the ancilla qubit in the appropriate basis can be classically processed and used for feedback operations on the system.the measurement of n-body observables such as multi-qubit stabilizer operators is an essential ingredient e.g. for error syndrome measurements in quantum error correction and quantum computing protocols [69,70,71].
The engineering of reservoir couplings and dissipative many-body processes enables novel directions for quantum state preparation [29], as recently also shown in an experiment with atomic ensembles [26].Combining dissipative time evolution with coherent Hamiltonian dynamics might allow one to explore novel physics such as nonequilibrium phase transitions in driven dissipative systems [72].In particular, the ability to implement e.g.master equations with multi-qubit quantum jump operators opens interesting perspectives for building quantum memories based on dissipation [33] or the demonstration of a novel form of quantum computing solely based on dissipation [32].
Star-type MS gate of the auxiliary ion with all system ions: It is also straightforward to realize an entangling gate between the ancilla ion and each system ion without creating pairwise entanglement between the system ions.This can be done by the sequence which is sketched in figure B1b.Here, the second inverse MS gate, which is applied after the single qubit flip of the auxiliary ion, cancels the initially generated entanglement between all pairs ions, which do not include the ancilla ion.
Sequence for a MS gate on 2 of n ions: The sequence for the implementation of the star-type entangling operation discussed in the previous paragraph can be used to realize a MS gate on two of n ions (see figure B1c).Such two-ion MS gate is an essential building block for the implementation of the two-qubit gate C(θ), which is needed for the dissipative simulation discussed in section 4. For instance, the sequence for a MS gate on only the auxiliary ion and the system ion #1 is given by U (0,1) MS (θ, φ) = U More involved decompositions for MS gates, where more than one or two ions are supposed to participate in or be excluded from the gate operation, can be constructed accordingly.In general, they will involve more "partial" MS gates and refocusing pulses, which might at some point render the alternative approach of hiding pulses on individual ions more suitable.

Appendix C. Gate Decompositions
Here, we provide decompositions of the two-qubit gates C i (θ), which are needed for the dissipative n-body dynamics as discussed in section 4, into MS gates and single-ion rotations.
For the simulation of n-body interactions with n = 4, 8, ... the gate C i (θ) of Eq. ( 15) can be decomposed as The two-qubit MS gate on the auxiliary ion and the i-th system ion U (0,i) MS (θ/2, π/2) can be realized via refocusing techniques -see eq.(B.3) in Appendix B.
It is straightforward to decompose the two-qubit gates C i (θ) for other values of n (as listed in table 2) accordingly.The nodes represent the n + 1 ions, lines between nodes i and j denote entanglement created between a pair {i, j} of ions during the application of MS gates.The short-hand notation ±MS/2 and ±MS/4 stands for U MS (±θ/2, φ) and U MS (±θ/4, φ), respectively.An operation j denotes a single ion pulse U (j) σ z (π) applied to the i-th ion; ions which have been exposed to such flip operations are labelled by a small circle, until they are flipped back into their original orientation.a) Gate decomposition for a MS on all ions except the n-th ion (cf.eq.(B.1)).b) Gate sequence for the creation of startype entanglement between the auxiliary ion (#0) and each of the n system ions.c) Gate sequence for a MS gate on two ions out of n + 1. Dashed lines correspond to entanglement, which is created in intermediate steps between an initially disentangled pair of ions {i, j}, if one (and only one) of the two ions (marked with a circle) has been previously flipped.

Figure 1 .
Figure1.Lattice spin models of interest for the gate-based ("digital") quantum simulation with trapped ions.(a) In Kitaev's toric code[35] spins located around vertices of a two-dimensional square lattice interact via four-body interactions ∼ σ x 1 σ x 2 σ x 3 σ x 4 , whereas spins around plaquettes experience z-type interactions, as e.g.∼ σ z 2 σ z 3 σ z 5 σ z 6 .(b) Small instance of a color code spin system, as proposed in[36].Here, spins are located on the sites of a three-colorable lattice interact via four-body plaquette interactions such as ∼ σ x 1 σ x 2 σ x 3 σ x 4 and ∼ σ z 1 σ z 2 σ z 3 σ z 4 .(c) Mesoscopic instances of spin models can be mapped onto linear chains of trapped ions, where the spin degree of freedom is encoded in (meta)stable electronic states.Coherent and dissipative time evolution can be simulated by sequences of highly parallel multi-ion Mølmer-Sørensen (MS) gates applied to all (or subsets of) ions, in combination with single-qubit rotations on individual ions and optical pumping of an ancilla ion.

Figure 3 .
Figure 3. (a) Graph representation of the two-body spin interaction Hamiltonian, which underlies the multi-ion MS gate(8).All pairs {i, j} of ions involved in the gate interact with equal strength (represented as links).(b) A (4 + 1) ion entangling MS gate applied to 4 system ions and an ancilla ion (index zero) can be used to coherently map the information about whether the 4 system ions are in a +1 (-1) eigenstate of the 4-body interaction term ∼ σ x 1 σ x 2 σ x 3 σ x 4 , onto the logical states |0 and |1 of the ancilla ion.In the Bloch sphere representation this mapping can be understood as a rotation of the ancilla qubit initially prepared in |0 around the x-axis.The rotation angle depends on the state of the system ions, and is chosen such that for any +1 eigenstate, as e.g.| + + + + , the ancilla qubit ends in |0 after the MS gate, whereas for -1 eigenstates such as e.g.| + + + − it is transferred to |1 .This mapping mechanism not only works for 4-body interactions, but can be used for general n.For n odd, the ancilla qubit is transferred to σ y eigenstates.

Figure 4 .
Figure 4. Generic gate decompositions for the simulation of coherent and dissipative dynamics via the Trotter expansion(7).(a) Coherent evolution according to a fourbody interaction Hamiltonian H α = E α σ x 1 σ x 2 σ x 3 σ x 4 for a time step τ is efficiently achieved in three steps: (i) First, an entangling MS gate U MS (π/2, 0) applied to the four system ions and the ancilla ion coherently maps the information on whether the system ions are in a +1 or -1 eigenstate of H α onto the logical states |0 and |1 of the ancilla qubit (cf.figure3).(ii) Second, a single-qubit gate exp(iφσ z 0 ) applied to the ancilla ion effectively imprints a phase −φ (φ) on all +1 (-1) eigenstates of H α .(iii) Finally, the initial mapping is reversed by an inverse MS gate, which disentangles the ancilla from the system ions, which have evolved according to exp(iφσ x 1 σ x 2 σ x 3 σ x 4 ).(b) Dissipative evolution, i.e. "cooling" of the system ions into the +1 eigenspace of H α mediated by Lindblad dynamics with four-body quantum jump operatorsc α = 1 2 σ z 1 (1 − σ x 1 σ x 2 σ x 3 σ x 4 ): Here, a two-qubit gate C(θ) = |0 0| 0 ⊗ 1 + |1 1| 0 ⊗ exp[iθσ yi ] is applied between the mappings (i) and (iii).In step (ii), all +1 eigenstates of H α are left invariant, whereas a spin flip σ y i with a θ-dependent amplitude is applied to one of the four system ions can convert -1 into +1 eigenstates.The angle θ controls the conversion probability and allows one to tune from small probabilities (θ 1, master equation limit) to unit pumping probability (θ = π/2).After the three steps, generally, the ancilla qubit is entangled with the system ions.(iv) Finally, the ancilla qubit is incoherently reset to its initial state |0 by optical pumping.This dissipative step enables to carry away entropy and "cool" the system qubits.

Figure 5 .
Figure 5. Cooling dynamics in a small-scale implementation of Kitaev's toric code.(a) Two-and three-qubit stabilizer operators for the boundary spins.(b) Two types of excitations in the model correspond to violations (eigenstates of eigenvalue -1) of x and z-type stabilizer operators.(c) Cooling dynamics due to pairwise "annihilation" of excitations or single excitations which are pushed out of the system.

Figure 7 .
Figure 7. Circuit for readout of the four-body stabilizer operator A = σ x 1 σ x 2 σ x 3 σ x 4 via a measurement of the ancilla qubit.The circuit for coherent simulation of a four-body spin interactions as discussed in section 3 (cf.eq.(13)) realizes the unitary operation U = exp(iπ/4σ z 0 ⊗ A).Thereby, the ancilla qubit initially prepared in |+ is coherently mapped onto the two σ y -eigenstates, depending on whether the four system qubits are in a +1 or -1 eigenstate of A. This information obtained from a subsequent measurement of the ancilla qubit in the appropriate basis can be classically processed and used for feedback operations on the system.

Figure B1 .
Figure B1.Gate sequences for the realization of entangling gates on subset of ions, by a combination of MS gates applied to all ions and refocusing pulses on individual ions.The nodes represent the n + 1 ions, lines between nodes i and j denote entanglement created between a pair {i, j} of ions during the application of MS gates.The short-hand notation ±MS/2 and ±MS/4 stands for U MS (±θ/2, φ) and U MS (±θ/4, φ), respectively.An operation j denotes a single ion pulse U