Engineering generalized Gibbs ensembles with trapped ions

The concept of generalized Gibbs ensembles (GGEs) has been introduced to describe steady states of integrable models. Recent advances show that GGEs can also be stabilized in nearly integrable quantum systems when driven by external fields and open. Here, we present a weakly dissipative dynamics that drives towards a steady-state GGE and is realistic to implement in systems of trapped ions. We outline the engineering of the desired dissipation by a combination of couplings which can be realized with ion-trap setups and discuss the experimental observables needed to detect a deviation from a thermal state. We present a novel mixed-species motional mode engineering technique in an array of micro-traps and demonstrate the possibility to use sympathetic cooling to construct many-body dissipators. Our work provides a blueprint for experimental observation of GGEs in open systems and opens a new avenue for quantum simulation of driven-dissipative quantum many-body problems.


I. INTRODUCTION
Providing a compact description of complicated manybody systems is a challenging task. Studies of equilibration of interacting quantum many-body systems have revealed that one can borrow statistical descriptions, taking into account conservation laws of equilibrating systems, to achieve that [1,2]. If we suddenly excite an ergodic system, for which energy is the only conservation law, steady state expectation values will be given by a Gibbs ensemble with temperature 1/β determined by the amount of energy that has been injected into the system with the excitation process [1]. Such a description by the single parameter -temperature -is, therefore, an incredible simplification for an interacting model with a priori exponentially many degrees of freedom. Similarly, generalized Gibbs ensembles (GGE) were proposed to describe local observables in steady states of systems with additional local conservation laws C i . The existence of additional conservation laws highly restricts the dynamics and prevents the system from thermalizing. GGEs have a form related to that of a Gibbs ensemble [3], but with additional Lagrange multipliers λ i associated with additional conserved quantities C i that all commute with the Hamiltonian, [C i , H 0 ] = 0. Exemplary models with macroscopically many local conservation laws, * These two authors contributed equally FIG. 1. Activating integrability. (a) Following an excitation of a non-integrable or nearly integrable system, the steady state exhibits thermal behavior with small expectation values of operators other than the Hamiltonian. (b) A nearly integrable system with a weak dissipative drive will show highly non-thermal behavior with a distribution p(En) approximately described by a generalized Gibbs ensemble (GGE) and large expectation values for conservation laws, e.g., C4 . (c) A one dimensional array of trapped-ions has interactions which decay as a power-law, and is thus nearly integrable for large exponents. We re-activate integrable effects by engineering the dissipation, which stabilizes a GGE as the steady state.
Applicability of GGEs was confirmed also experimentally in a cold atoms setup [26] where, up to some time, a closed and integrable system can be prepared. Ref. [26] showed that GGEs for a Lieb-Liniger model can provide an accurate description of an interacting trapped 1D Bose gas. However, it is very difficult to simulate integrable systems due to their fine-tuned nature: adding practically any other terms to an integrable model will break integrability and cause eventual thermalization [27][28][29][30][31][32][33][34][35][36]. Therefore it was believed that in realistic systems traces of integrability can be seen only in the transient dynamics [37][38][39][40][41][42][43][44][45][46] while the steady state is always thermal due to realistic integrability breaking terms.
In recent works [47][48][49], Lange and Lenarčič demonstrated that properties related to integrability are not as fragile as previously believed: if one weakly drives an only approximately integrable system and at the same time allows it to cool via a weak coupling to the environment, the system will nonetheless relax to a steady state approximated with a generalized Gibbs ensemble, supplemented with a small correction δρ, The major difference to the closed strictly integrable setup, Eq. (3), is that here the Lagrange multipliers λ i are determined by the integrability breaking perturbations themselves, through a stationarity condition [47][48][49] ∂ t C k ≈ tr C k L p e − i λiCi tr e − i λiCi ! = 0, ∀k.
Here L p denotes the Liouville operator corresponding to perturbations which weakly break the integrability, and as a consequence the conservation laws, while driving and cooling the system. Such a setup is much more versatile because it does not require the fine-tuned perfect integrability and at the same time allows for the engineering of GGEs through a particular choice of perturbations. A remarkable consequence is that one can stabilize steady states with large expectation values of nearly conserved operators C i , if a large corresponding λ i is established. For example, in solid-state spin chain materials, approximately described by an XXZ model coupled to phonons, a weak laser driving could stabilize steady states with huge heat and spin currents, since these are (partial) conservation laws of the XXZ model [47]. Alternatively, driving and openness can be provided by Markovian dissipative processes [48]. As shown in Fig. 1, the latter could be experimentally realized with trapped-ion platforms where long-range interaction between ions inevitably breaks integrability. By adding engineered dissipative processes, a GGE ensemble would be stabilized. As a consequence, the distribution p(E n ) = exp(− i λ i n|C i |n )/tr e − i λiCi over the eigenstates of H 0 |n = E n |n would differ from a thermal one p(E n ) = exp(−βE n )/tr e −βH0 . While p(E n ) is essentially impossible to measure experimentally, the non-thermal nature of ρ GGE is more easily detected through possible exceedingly large expectation values of conserved operators C i of parent XY or transverse-Ising Hamiltonian.
In this work, we address the latter option and discuss an implementation based on ion-trap technology [50][51][52]. Controllable coherent couplings in ion-trap systems have been widely used for the realization of spin models [53,54], and numerous milestone experiments have been conducted on a variety of systems [55][56][57][58][59][60][61][62][63]. State-of-theart Paul traps [51,59,62], Penning traps [61,64,65], and micro-traps [66][67][68][69][70] offer a rich toolbox of couplings suitable to engineer coherent Hamiltonian, as well as dissipative interactions. Sympathetic cooling based on mixedspecies ion chains is well-studied in Paul traps [71][72][73] where it is used to remove entropy from the motional modes of the ions. Here, we realize a driven-dissipative dynamics consisting of a spin Hamiltonian, in combination with one-and two-body dissipation. The dissipators are engineered combining tunable carrier and sideband couplings with repumper drives or sympathetic cooling as sources of dissipation. To tightly confine the motional modes to the interacting particles, we propose a novel mixed-species mode engineering technique that can be realized in micro-trap arrays. The resulting dynamics stabilizes a steady state approximately described by a generalized Gibbs ensemble, despite different integrability breaking terms.
The paper is organized the following way: In Sec. II we introduce one choice of a Hamiltonian and Lindblad operators that could be realized in a trapped ion experiment. In Sec. III we present numerical results and discuss experimental signatures and means to measure that a GGE approximates the stabilized steady state. In Section IV-V we present the engineering of the elementary dissipators. We then scale up the interactions in Sec. VI and Sec. VII, where we discuss the mode engineering in arrays of micro-traps.

II. MODEL
The theory of activating integrability and engineering steady states described by generalized Gibbs ensembles in realistic systems is generic and applies to different systems approximately described by an integrable model, e.g., transverse field Ising or XXZ Heisenberg chain, Lieb-Liniger or Tonks-Girardeau Bose 1D gas. We choose an integrable model that is closest to the state-of-the art trapped ions setups. We consider the XY-Hamiltonian in the presence of a magnetic field h, which belongs to the class of non-interacting integrable models. Such a Hamiltonian can be implemented by using standard techniques developed in the trapped-ion field. In contrast to Ref. [59], we rotate the spin axes by π/2 around the y-axis. The resulting YZ-model will allow us to facilitate an experimental implementation of Lindblad terms. We consider An alternative realization based on the XY-Hamiltonian in combination with sympathetic cooling is presented in Sec. V. In traditional setups with trapped ions, the coupling between spins that are d sites apart actually decays as j+d with α ∈ [2,3]. This is one inevitable source of integrability breaking since such Hamiltonian can no longer be diagonalized via a Jordan-Wigner transformation, neither is it Bethe-Ansatz solvable. However, if the decay is fast enough one can consider such a system as nearly integrable. In our analysis we will take into consideration only the leading contribution H 1 alone would thermalize the system, however, nonthermal steady states approximated by GGE can be achieved when a weak coupling to Lindblad nonequilibrium baths is added. Here, we will consider the homogeneous bulk dissipators of two types, a = 1, 2, with Lindblad operators at site j Here is a projection on the spin-down state at site j.
The steady state density matrix ρ ∞ is determined bẏ where L 0 is a dominant term in the Liouvillian, while L p = L u + L m captures unitary and Markovian perturbations, Despite the fact that the underlying model H 0 in Eq. (6) is non-interacting, the next-nearest neighbor interaction H 1 and our choice of Lindblad operators hinders analytical solvability and requires a numerical solution. Since our Lindblad operators have local nature, we also cannot use recently proposed hydrodynamic description [74]. For a review of the theory of weakly driven nearly integrable systems, developed in [47][48][49] see App. A. In the next section we go on to establish that this model does stabilize a steady state approximately described by a GGE ensemble, despite different sources of integrability breaking.

III. EXPERIMENTAL SIGNATURES
The experiment we are proposing would aim to show that a highly non-thermal steady state, described with a generalized Gibbs ensemble, can be stabilized in a nearly integrable model given by H 0 + H 1 , Eqs. (7,8), if weakly driven ( 1) with dissipation Eqs. (10,11). To show that a GGE is stabilized by driving, we compare in Fig. 2 the steady state expectation values obtained with the exact density matrix, Eq. (12) or with a GGE, Eq. (14). A qualitative agreement confirms that a GGE is stabilized. The first notable conclusion from the numerical analysis is that in the limit of small driving, 1, a truncated generalized Gibbs ensemble parametrized with only N c = 4 Lagrange parameters can capture expectation values of local observables we are interested in, while other more complicated conservation laws can be neglected. In comparison, a full density matrix requires 4 N parameters in the system size N . This shows that a description in terms of truncated GGEs is extremely compact. For a more thorough comparison see App. A. Another interesting observation is that the nextnearest neighbors coupling H 1 , Eq. (8), does not have a strong impact on the steady state. While in a closed setup H 1 is crucial as it dictates relaxation towards a thermal state, in our setup it is dominated by Lindblad terms, see App. A for details. Therefore, we neglect it in results presented in the main text.
A physically most interesting consequence of stabilizing a steady state approximately described with a GGE is that the expectation values of conservation laws, can be much larger than in a thermal state, tr C k e −β0 / tr e −βH0 . In this respect conservation laws are measurably distinct from other operators: a generic observable O will have a much smaller expectation value H0 and C4 as a function of relative driving strength γ, Eqs. (10,11). A qualitative agreement between the exact ρ∞ and ρtGGE confirms that a generalized Gibbs ensemble is stabilized. A large C4 , which could not be obtained from a nearly infinite temperature steady state consistent with a small H0 , indicates that a highly non-thermal GGE is stabilized. Parameters: Jy = h = 1, Jz = 0.1; N = 10, Nc = 4 for ρtGGE and N = 6, = 0.01 for ρ∞. than a conservation law, tr[Oρ GGE ] tr[C k ρ GGE ], given they have the same norm tr O 2 = tr C 2 k . The conservation law C k will show a particularly large expectation value if driving is such that it stabilizes a GGE with a large corresponding Lagrange parameter λ k . Which λ k = 0 depends on the symmetry of the driving, i.e., the Lindblad operators [47,48,75]. A strong response of conservation laws to weak driving could have practical implications, such as heat and spin pumping [47] in spin chain materials, but can also serve to detect that a GGE has been stabilized despite the integrability breaking terms. In the following we discuss ways to detect large expectation values of conserved (or partially conserved) operators, which could not be possible in a thermal state.
One possibility is to compare expectation values of observables which do or do not overlap with conservation laws. If driving stabilizes ρ GGE with a large Langrange multiplier λ i associated with the conservation law C i , observable O that has a nonzero overlap with C i , tr[OC i ] = 0, will show a large expectation value. Let us, for example, consider σ y j−1 σ x j σ y j+1 ≡ YXY and σ y j−1 σ y j σ x j+1 ≡ YYX . At least at small J z /J y 1, where H 0 ≈ J y σ y j σ y j+1 +hσ x j , it is easy to estimate thermal expectation values XYX th , XXY th via expansion ρ th = e −βH0 /Z ≈ (1−βH 0 +β 2 H 2 0 /2+. . . )/Z: a nonzero YYX th is dominantly coming from the 2nd order in β, while YXY th from the 3rd order. In a thermal state one would therefore expect YXY th YYX th . However, our numerical result in Fig. 3 shows YXY YYX . This observation is a clear sign that the steady state is not thermal. The large expectation value of YXY is a direct consequence of the fact that this operator is part of C 4 , Eq.(16), therefore its expectation value has also a linear contribution in λ 4 in the expansion, since driving stabilizes a GGE with λ 4 = 0. In our setup, besides the Hamiltonian H 0 (H 0 = C 2 according to the notation we use, see App. A), the next simplest local extensive conservation law C i = j c (i) j that shows a strong response to the dissipative driving is wherez = y andȳ = z. Fig. 2 shows the comparison of H 0 = C 2 and C 4 as a function of the relative strength γ, Eqs. (10,11), obtained from a GGE calculation on N = 10 sites or a full steady state density matrix on N = 6 sites. While H 0 (γ) shows rather small values, C 4 (γ) is bigger. That would never be the case, if C 4 (γ) was evaluated with respect to (an almost infinite temperature) thermal ensemble, which reproduces small H 0 (γ). This observation alone suggests that a non-thermal steady state is stabilized. In order to quantify how non-thermal the steady state is and to select optimal parameters, we introduce the calculated with respect to the exact steady state, ρ x = ρ ∞ , or with the truncated GGE, ρ x = ρ tGGE . For calculations with ρ x = ρ ∞ we define ρ th as a thermal state with respect to H 0 , Eq. (6), with temperature determined from the condition tr For calculations based on ρ x = ρ tGGE the temperature in ρ th is calculated using a Gibbs ensemble ansatz with H 0 as the only conservation law.
In the following we focus on the operator O = C 4 . In the absence of Lindblad driving D 1 (2) , the ratio equals η C4 = 0 due to integrability breaking power-law decay of interactions in the Hamiltonian. In the presence of a weak Lindblad drive, on the other hand, the steady state 0.

4.
8. can be highly non-thermal. Figure 4 shows η C4 as a function of anisotropy J z /J y , obtained from the exact density matrix ρ ∞ on N = 6 sites at = 0.01 and from a ρ tGGE on N = 10. The dependence on J z /J y suggests that the experiment observing a highly non-thermal steady state (η C4 1) achieved by a weak driving should operate at a small J z /J y or even at J z = 0 (corresponding to the transverse-field Ising case).
Why η C4 is so large for J z = 0 (and small J z ) can be reasoned by looking directly at the expectation values of H 0 and C 4 , Fig. 5. H 0 ≈ 0 is almost zero, suggesting an almost infinite temperature state (β ≈ 0), which would imply C 4 ≈ 0. On the other hand, the observed C 4 ∼ 1 is actually large. This observation, together with a good agreement between ρ ∞ and ρ tGGE , clearly shows that a GGE with a large λ 4 is stabilized, despite  β ≈ 0. Measuring H 0 and C 4 would thus serve as strong affirmation of our theory.
In Fig. 6, we show that at mild anisotropy, for example, J z /J y = 0.9, also a magnetic field h helps to prepare a more non-thermal state.

IV. DISSIPATION ENGINEERING
Having shown that our dissipative driving can stabilize a steady state described by a GGE, we now discuss the implementation of the desired dynamics in a trapped-ion setup. To engineer suitable dissipative interactions, we combine coherent couplings with sources of dissipation such as induced spontaneous emission and sympathetic cooling. We use these tools to engineer the desired oneand two-body jump operators and verify their action numerically.
In this section, we assume that the YZ-Hamiltonian in Eq. (7) is implemented. This allows us to engineer the dissipation in Eqs. (10)-(11) between the levels |↑ , |↓ . As these are eigenstates of (stimulated) spontaneous emission, it is possible to use repumper beams as sources of dissipation. In Sec. V, in turn, we assume that the XY-model in Eq. (6) is realized. This requires the dissipation to be engineered between the eigenstates of σ x , as is demonstrated using sympathetic cooling.

A. Setup
To implement the one-and two-body jump operators in Eqs. (10) and (11), we consider a system of trapped ions coupled through motional modes. To simplify the discussion, we start by considering a minimal instance consisting of two ions indexed 1 and 2, and a motional mode a, and generalize to more ions in Sec. VI-VII. As is shown in Fig. 7, each of the ions is assumed to have two Setup for engineered single-body and two-body dissipation. We consider two trapped ions, j = 1 and j = 2, with two stable ground levels, |↓ and |↑ , and three excited levels, |e , |r , and |q . The ions are subject to coherent couplings (solid lines) and dissipation (dotted lines). (a) A weak drive excites both ions from |↓ to |e (strength Ω, detuning ∆). The transition from |e to |↓ is coupled to the motional mode a (phonon detuning δ) by a red-sideband interaction (coupling constant g). The second ket denotes motional excitation and is dropped when in vacuum. (b) We realize single-body decay from |↑ to |↓ and from |e to |↑ on ion 1. This is done by optical pumping using tunable repumper beams (Rabi rates Ω rep,↑ and Ωrep,e), via the unstable levels |r and |q (decay rates Γ ↓r and Γ e↑ ).
stable ground levels, |↑ and |↓ , and three excited levels, |e , |r , and |q . The motional mode a is assumed to be cooled to the ground state, |0 . The free Hamiltonian of this system is given by Here we introduce a phonon detuning δ and an ionic detuning ∆, assuming that we work in a suitable rotating frame with respect to the fields to be introduced below. We will use level |r to realize the single-body decay in Eq. (10) and levels |e and |q in combination with mode a for the two-body dissipation in Eq. (11). Level |q will be used to add a dissipation channel from level |e to level |↑ . The ions are excited from |↓ to |e using a weak "carrier" drive with a Rabi frequency Ω. In addition, levels |↑ and |e of ion 1 are excited to |r and |q by coherent "repumper" beams, H rep,e = Ω rep,e 2 |q 1 e| + H.c.
with Rabi rates Ω rep,↑/e . The coupling between ions 1 and 2, needed to engineer two-body dissipation, is medi-ated by a common motional mode, with creation (annihilation) operator a † (a). This phonon mode is coupled to the transition from |e to |↓ by the red-sideband interaction with a coupling constant g.
In addition to the above coherent interactions, the excited level |r is assumed to be inherently unstable and to decay to level |↓ by spontaneous emission, which can be described using the jump operators, The excited level |q , in turn, is assumed to decay to |↑ , To describe the joint dynamics of the ions and the phonons, we use the following notation: the state of the system is described by two kets, where the first ket denotes the internal state of the ions, e.g., |↓↓ = |↓ 1 |↓ 2 . Motional excitations are denoted by a second ket, e.g., |↓↓ |1 , which is dropped when the motion is in the ground state.

B. Single-body dissipation
To realize single-body dissipation, we employ standard optical pumping, combining excitation from |↑ to |r by H rep,↑ , Eq. (20), and decay from |r to |↓ by spontaneous emission L ↓r,j . The effective jump operator [76] for the decay of level |↑ to |↓ through |r is thus, after elimination of level |r , given by We thereby realize the desired jump operator L . Using individual addressing techniques, this process can be made site-specific. The decay rate γ 1 can be tuned by varying Ω rep,↑ , assuming it to be much smaller than the natural linewidth of level |r , Γ ↓r Ω rep,↑ . Note that, while here we have only assumed the desired dissipation channel from |r to |↓ , additional decay processes from these levels could be described by the same method.
To engineer the two-body dissipation in Eq. (11) in Sec. IV C below, we will also rely on an induced spontaneous emission process from |e to |↑ . Here we assume that we can realize the YZ-Hamiltonian in Eq. (7). Alternatively, in the presence of the XY-Hamiltonian in Eq. (6), sympathetic cooling can be used as a source of dissipation, as is described in Sec. V. (|e ↓ + |↓ e ) |0 by the drive Ω. |ψe |0 is strongly coupled to the motion-excited states |↓↓ |1 by the sideband coupling g, which is enhanced by a factor √ 2 due to constructive interference. For ∆δ = 2g 2 , the lower dressed state of |ψe |0 and |↓↓ |1 (indicated in blue) is in resonance with the drive and hence rapidly excited from |↓↓ |0 . Through its contribution from |ψe |0 , it decays to |↑↓ |0 by spontaneous emission Γ. These resonant couplings form an effective decay process from |↓↓ |0 to |↑↓ |0 , at an enhanced rate γ2. (b) Undesired process. Also state |↑↓ |0 is excited by the drive, to an ion-excited state |e ↑ |0 . The sideband coupling couples |e ↑ |0 to |↑↓ |1 at a coupling constant g. Their dressed states are thus out of resonance and only weakly populated by the drive.
To realize optical pumping from |e to |↑ , we couple |e of ion 1 to the unstable level |q (total decay rate Γ q ) using the repumper in Eq. (21). Together with the decay from |q to |↑ in Eq. (24), this realizes Again, the decay rate Γ is tunable through the strength of the corresponding repumper beam Ω rep,e .

C. Two-body dissipation
We now turn to the two-body dissipation in Eq. (11). Compared to the single-body dissipation in the previous section, the operator L (2) j is more complicated to engineer, but itself sufficient to realize a highly non-thermal GGE (see Fig. 2 at γ = 1). For our minimal instance of two ions, the operator reads L The action of this operator can be understood as a raising on spin 1, S + 1 = |↑ 1 ↓|, conditioned on the state of spin 2.
We will now engineer the desired two-body dissipation based on the assumptions of weak driving, Ω 2 {Γ 2 , g 2 }, and strong coupling, Γ 2 g 2 . In this regime, the ground state |↓↓ |0 is weakly excited by H drive to the excited state |ψ e |0 = 1 √ 2 (|e ↓ + |↓ e ) |0 , which comprises of a superposition of excitations of both ions. Further excitation to the double-excited state |ee can be neglected, as we will see further down. We engineer the desired mechanism using the couplings of |ψ e |0 : The ion-excited state |ψ e |0 is coupled to the motionexcited state |↓↓ |1 by the red-sideband coupling H int . Due to constructive interference between the excitation of the two ions, the corresponding rate is given by √ 2g. The Hamiltonian for the coupled subspace is and illustrated in Fig. 8 a). Based on the assumption of a weak drive compared to the rapid dynamics of the excited subspace, we make a separation of timescales and first regard H e,↓↓ alone, without the drive: Provided strong coupling, the excited states of the strongly coupled subspace hybridize and form dressed states |ψ ± at detunings Setting the ionic and the motional detunings to ∆δ = 2g 2 (e.g., ∆ = δ = √ 2g) brings the lower dressed state in resonance with the drive Ω, i.e., ∆ − = 0. As a consequence, |↓↓ is resonantly excited to |ψ e which in turn decays to |↑↓ at a rate Γ. This results in the required effective decay from |↓↓ to |↑↓ , mediated by the resonant lower dressed state, |ψ − . In App. B we present a full microscopic derivation that yields the two-body Lindblad operator, Eq. (11), Using second-order perturbation theory in Ω/Γ the effective decay rate is given by, The decay rate γ 2 , and hence the relative strength of the single-and two-body dissipation, can thus be adjusted by varying Ω and Γ(Ω rep,e ). γ 2 is ultimately limited by the linewidth Γ q of level |q , which is involved in constructing the engineered single-body decay in Eq. (26). Compared to γ 2 , effective decay processes mediated by |ee can be neglected: |ee forms a coupled two-excitation subspace with |ψ e |1 and |↓↓ |2 (couplings √ 2g and 2g). While for the above parameter choice, the drive from |ψ e is in resonance with the two-excitation dressed states, the coupling rates of the so mediated effective process only enter to fourth order in perturbation theory and are thus negligibly small. Also AC Stark shifts arising from the weak off-resonant excitation of |↑↓ (cf. App. C) can be safely ignored in the considered parameter regime.
Another imperfection inherent to the scheme is given by the population of the excited level |↑ e which is offresonantly excited from |↑↓ by the drive Ω, as illustrated in Fig. 8 b). For perfect individual addressing of the first ion by H rep,e , |↑ e is steadily populated. Using adiabatic elimination, it is possible to estimate the steady-state population of an excited state such as |↑ e [77,78], which scales Population of the excited state |↑ e is thus expected to be largely suppressed, as we numerically confirm in the following.
We verify the engineering of the two-body dissipation and assess its performance numerically by simulating the dynamics given by the master equatioṅ We use this to obtain i) the optimal parameter choice, ii) the extent of unwanted population of the excited level |↑ e , and iii) the timescales of the dissipative compared to the unitary dynamics. We assume that the system starts from |↓↓ |0 and optimize the fidelity of the state |↑↓ |0 after a chosen time, t opt = {50, 100, 200}/g by the choice of the parameters Γ and Ω. The result is plotted in Fig.  9. From the initial state |↓↓ , the system evolves to fidelities with |↑↓ of Generalization of the couplings to the x-basis. (a) To engineer two-body dissipation in the x-basis, we use couplings between the states |± = (|↓ ± |↑ )/ √ 2, and |e . This is achieved by coupling the levels |↓ and |↑ coherently to |e . (b) Dissipation in the x-basis is facilitated by sympathetic cooling of the motion. The single-body decay process from |+ to |− is engineered by excitation from |+ by a weak drive (Ωrep,+) to an auxiliary level |r , a sideband coupling of the transition |r → |− to a motional mode b (g b ), and sympathetic cooling of b (κ b ). Decay from |e to |+ (not shown) is engineered accordingly using a motional mode c, as is described in the text. exponential form that is described by the effective jump operator in Eq. (29). This represents a close approximation to the desired dynamics.
The optimal parameters fulfill the conditions for weak driving and strong coupling, Ω 2 Γ 2 g 2 , which are the assumptions used in the above analysis. The residual population in |↑ e is found to be P ↑e ≈ {0.02, 0.008, 0.003} and thus indeed negligible.
The achievable values for γ 2 should be compared with the typical coupling strength in realizations of the spin models of about J x/y /(2π) ∼ 100 Hz [59,60], and the corresponding timescale of τ spin ∼ 1 ms. From the above numbers it can be seen that the effective dissipation rate can be tuned to values within the same order of magnitude as the coupling constants of the spin model. Obtaining smaller values for γ 2 -and thus making the dissipation into a perturbation as assumed in Sec. II-III -is in turn achieved by choosing weaker repump and driving rates, Γ and Ω.

V. GENERALIZATION OF THE COUPLINGS
TO THE x-BASIS Spin models along arbitrary directions, with and without anisotropy, can be realized on trapped-ion platforms, such as Paul traps and Penning traps, and also in micro-traps [53]. The majority of the available setups, however, support XY-Hamiltonians without anisotropy [59]. In the preceding section, we assumed the less common YZ spin Hamiltonian which enabled us to engineer the dissipation in the z-basis, the eigenbasis of decay by spontaneous emission. As an alternative to such implementation, we can utilize the more standard XY-Hamiltonian in Eq. (6). We should note, however, that in order to stabilize a non-trivial (ρ ∞ = 1) and distinguishingly nonthermal steady state, an engineered form of Lindblad operators with proper symmetries must be used. For example, the XY-Hamiltonian in combination with the dissipation in the z-basis would not result in a large C 4 . This problem is resolved by using the XY-Hamiltonian in combination with dissipation in the rotated x-basis, the engineering of which is presented below.
We point out that a simpler alternative is given by using the XY-Hamiltonian in combination with any type of Lindblad dissipation. Such setting would yield a nontrivial dynamics towards a possibly trivial steady state, where time evolution is approximately captured with a time-dependent GGE [48]. However, the non-thermal features might not be very pronounced in a generic setup, therefore we focus here on a setup with non-trivial steady states with clearly non-thermal nature.
We now engineer the dissipation in the x-direction which can stabilize a steady-state GGE, when combined with the XY Hamiltonian. Formally speaking, the change from the z-to the x-basis amounts to replacing the eigenstates of σ z , {|↓ , |↑ } with the eigenstates of σ x , |± = (|↓ ± |↑ )/ √ 2, i.e., |↓ → |− , |↑ → |+ , in all steps of our previous derivation performed for the dissipation in the z-direction. Physically, using |± is realized by coupling to both transitions |↓ ↔ |e and |↑ ↔ |e coherently. The resulting interactions, illustrated in Fig. 10 a), read On the other hand, decay by spontaneous emission, as utilized in the previous sections naturally occurs in the z-basis, {|↓ , |↑ }. In contrast to this, we now need to engineer sources of dissipation in the x-basis, which replaces Eqs. (23), (26) in the z-basis. In the following, we demonstrate how to achieve this using decay of excitations via the motional degree of freedom by sympathetic cooling. As shown in Fig. 10 b), to implement the single-body decay in Eq. (33), we excite |+ to the auxiliary level |r by a repumper H rep,+ = Ω rep,+ |r 1 +| + H.c.
The excitation to level |r is transferred coherently to an auxiliary motional mode b using a sideband interaction, with a coupling constant g b . Mode b is subject to sympathetic cooling which realizes the jump operator Adiabatic elimination of b leads, for g b κ, to the desired decay channel with a rate Γ −r = g 2 b /κ, in analogy to Eq. (23). We realize the decay L +e = Γ +e |+ 1 e| (cf. Eq. (26)) in a similar fashion, utilizing a motional mode c, which is subject to sympathetic cooling. Here, we couple |e |0 c to |+ |1 c by a sideband drive (g c ), which then decays to |+ |0 c by sympathetic cooling (κ c ), resulting in an effective decay rate Γ +e = g 2 c /κ c . Involving level |q is not necessary. Carrying out the same analysis as in Sec. B, we obtain the effective operator with a tunable decay rate γ 2,x = 4Ω 2 /Γ 2 +e . We have thus realized the desired two-body dissipation in the x-basis in Eq. (34) by means of sympathetic cooling.

VI. SCALABLE IMPLEMENTATION
Next, we discuss how to scale the mechanisms discussed in Sec. IV-V to larger numbers of ions. For a scalable implementation of our scheme, we assume a chain of N ions (with even N ) and a level structure similar to Sec. IV A. The physical system for the scalable implementation of two-body dissipation in Eq. (11) is shown in Fig. 11. The implementation of such system based on ion micro-traps is presented in Sec. VII. Here we seek to implement interactions on all pairs of ions, such as {2j − 1, 2j} and {2j, 2j + 1}, as illustrated in Fig. 11 a). However, care has to be taken to avoid interference effects of the coherent couplings in the overlapping region, i.e., here ion 2j. We achieve this by devising two independent coupling configurations to mediate the engineered decay on the two different groups of ions, {2j − 1, 2j} and {2j, 2j + 1}, as can be seen from Fig. 11 b)-c). We assume each ion to have two (meta-) stable excited levels, |e and |f , which are selectively addressable using, e.g., polarization selection rules.
To implement two-body dissipation, we again use tunable optical pumping of |e and |f to |↑ , in analogy to Sec. IV C, Eq. (26). Utilizing different individually addressed repumper beams for "odd" ions 2j −1 and "even" ions 2j, we realize Odd ions 2j − 1 thus decay from |e to |↑ , whereas even ions 2j decay from |f to |↑ , both at an equal rate Γ.
For the continuous and measurement-free interrogation of the system, we use two sets of coherent drives, H drive,e = Ω 2 N/2 j=1 |e 2j−1 ↓| + |e 2j ↓| + H.c., (45) H drive,f = Ω 2 N/2 j=1 |f 2j ↓| + |f 2j+1 ↓| + H.c., (46) coupling the ground level |↓ to the excited level |e or |f , as well as sideband interactions, These realize coupling configurations, by which the transition |e ↔ |↓ (|f ↔ |↓ ) of any pair of ions {2j −1, 2j} ({2j − 1, 2j}) is coupled to a localized motional mode a 2j−1,2j (a 2j,2j−1 ). As a result, following the recipe in Sec. IV C, we realize jump operators acting on pairs of ions over the whole chain, In the second step, these operators are brought back into the form of Eq. (11). Making the association γ 2 = γ, we have thereby engineered the desired two-body dissipation in a scalable manner. Single-body dissipation in Eq. (10), is again realized -now for the whole chain -following the recipe in Sec. IV B: Using locally addressed repumper beams to an unstable level |r for each individual ion, we achieve local jump operators Associating γ 1 = (1 − γ), we have thus realized the desired single-body dissipation (Eq. (10)) for all ions in the chain.
The implementation of generalized dissipation in the x-basis, such as those in Eqs. (33)- (34) in Sec. V can be scaled up in an analogous manner.

VII. NORMAL MODE ENGINEERING IN AN ARRAY OF MICRO TRAPS
In the following, we discuss the physical implementation of the scalable setting detailed in Sec. VI based on ion micro-traps [66][67][68][69][70]. To implement the desired mode structure, with localized modes subject to dissipation, we employ a novel approach to mixed-species normal mode engineering.
In an array of micro-traps individual control over the electric potential at each trap site allows us to engineer normal mode spectra suitable for implementing the desired operators. Previously the use of a special arrangement of the transverse frequencies of individual traps along a string has been proposed for encoding two-body bosonic gauge fields [68]. Here we consider similar ideas to construct phonon modes localized to triplets of neighboring ions (that is any ion together with both of its nearest-neighbors), and with the use of mixed species of ions effectively encode two-body as well as one-body operators. The use of two species of ions allows us to achieve this using mode engineering along only one axis of vibration instead of two transverse axes.
The calculation of the normal modes for a system of ions in an array of micro-traps [70,79] involves first deter-mining the equilibrium positions in the combined electric potential due to the trap electrodes and the Coulomb repulsion from the other ions. A Taylor series expansion about these positions up to second order then yields the effective harmonic potential experienced by each ion. The resulting Hamiltonian is then diagonalized to extract the eigenvalues, which give the frequencies of oscillation, and the eigenvectors, which give the relative amplitude of oscillation of each ion in any given normal mode. Since the trapping mechanism in Paul traps and Penning traps differs only in the radial directions, we focus here on the axial modes since the treatment would then be independent of the kind of trap for the discussion below.
To illustrate the idea we first consider a chain of a single species of trapped ions, each with mass m. The spatial separation between neighboring traps is given by d and the potentials are designed so that the trap frequencies ω k = ω 0 + (k − 1)∆ T increase linearly with the trap site k. If the frequency difference ∆ T is much larger than the two-ion exchange frequency, defined by Ω ex = e 2 / 4π 0 mω 0 d 3 , we get for this system of coupled ions a spectrum of normal modes close to that of the non-interacting system. That is, the result is a set of modes with oscillation frequency near ω k and participation mostly from the ion at that trap k (as well as ions at the nearest-neighboring traps, k − 1 and k + 1). Now consider a chain of alternating ionic species of nearly equal mass in these traps. One of these species serves as a coolant ion while the other one is used to encode the spins through two suitable internal states, and is from here on called the 'spin ion'. The condition for nearly equal mass is useful for efficient sympathetic cooling of the spin species through the coolant species. As discussed above, the resulting normal mode structure consists of modes effectively localized to triplets of ions (except, of course, at the edges) but now there exist two kinds of modes -one where the central ion is the coolant ion with two neighboring spin ions, and the other where the central ion is a spin ion between two coolant ions. Fig. 12 a) shows this for a chain of alternating 40 Ca + and 43 Ca + ions arranged along the axial direction of the micro-trap array. The frequencies for the uncoupled arrangement of traps and for the coupled system are shown in figure 12 b). Here ω 0 = 2π × 300 kHz, ∆ T = 2π × 37.5 kHz, and d = 30 µm. Since the participation of the ions drops roughly exponentially with the distance from the central trap we can assume almost no participation from other ions of the same species (since these lie two sites distant). This behavior can be seen from Fig. 12 c), where the amplitudes of motion of each ion in the axial modes is plotted. Note that each mode is assigned a color so that it is clear which ions dominate the oscillation in a particular mode, and at what frequency. Modes a which are localized to two spin-ions (and one coolant ion) allow us to engineer the two-body operators while modes b which are localized to one spinion (and two coolant ions) allow for implementing onebody operators. To match the notation of the Sec. VI, we associate the even-numbered trap sites with the spin ions by 2k → j (e.g. a 2,4 → a 1,2 ), thereby excluding the odd-numbered trap sites containing the coolant ions. Employing sympathetic cooling of these modes, dissipation can be engineered in the σ x -eigenbasis, following the recipes in Sec. V.
Note that the laser-ion couplings in Sec. VI, Eqs. (44)-(49) require addressing of pairs of ions coherently with the same detuning ∆. As discussed in Sec. IV C, ∆ needs, however, to be matched to the phonon detuning δ, which differs from ion to ion because to the trap frequency offset. The resulting mismatch can be compensated by local AC Stark shifts on |e and |f , which can be generated by individually addressed lasers.

A. Alternative realization in linear ion traps
In conventional bulk ion traps, the desired dynamics can be implemented in a stepwise manner. Here we consider a sequential realization based on delocalized motional modes in combination with local addressing techniques [59,[80][81][82][83]. We assume all ions to be detuned by a constant amount with respect to the coupling configurations in Sec. VI. For N ions, we now consider N timesteps. During each step, we direct a pair of individual addressing lasers on a pair of ions. This beam shifts the transitions of the ions into resonance with the carrier and sideband couplings in Eqs. (44)- (49). We thereby pairwise realize the desired dynamics on ions j and j + 1, leaving the other ions uncoupled. In the next step, the individual addressing laser is shone onto another pair of ions, j + 1 and j + 2, and so forth. Provided short modulation times for the lasers, the timeframes for the ions may be reduced to stroboscopic length, resulting in a "Trotterized" realization of the desired dynamics.

VIII. CONCLUSION AND OUTLOOK
We have presented a scheme suitable to revive the effects of integrability in a controllably driven and open setup, where the underlying Hamiltonian dynamics is only approximately integrable. The scheme is based on weak couplings to Markovian baths in combination with nearly integrable quantum spin Hamiltonians, ingredients which are readily available in state-of-the-art trapped-ion setups. In addition, we present a novel technique for the engineering of motional modes in an array of mixed-species micro-traps, which support the realization of the desired dynamics.
Our numerical analysis shows that despite different sources of integrability breaking due to long-range interactions in the Hamiltonian and openness itself, a steady state is realized that cannot be modeled as a thermal ensemble. Instead, approximate expectation values of local observables can be obtained from a generalized Gibbs ensemble and we identify the experimental signatures which reveal that. We presented results for a rotated XY and transverse-field Ising Hamiltonian; however, the same Lindblad operators would activate a GGE in the interacting XXZ Heisenberg chain.
While our goal has been to engineer Lindblad dissipators which stabilize a non-trivial and highly non-thermal steady state by weak driving, a simpler task would be to consider a non-trivial dynamics towards a trivial (thermal, maximally mixed or empty) steady state. As has been shown in Ref. [48], the dynamics of a (nearly) integrable system weakly coupled to arbitrary Lindblad baths can be approximately described with a time-dependent GGE. An example of this is atom loss in cold-atom setups, whose effects have been observed [84] and theoretically addressed [85] very recently. The simplest combination of the strategy presented above would be to follow the dynamics of an XY-Hamiltonian, Eq. (6), in the presence of single-body decay, Eq. (10), alone.
Beyond this work, our dissipation engineering strategies open the door to experiments that will shed light on novel phenomena in open quantum systems. For example, the dissipators presented here have been used to study many-body localization through the perspective of open systems [86]. Our novel mixed-species mode engineering techniques based on arrays of ion micro-traps hold promise to become a powerful tool in quantum simulation. Here, sympathetic cooling is not only useful to reduce the entropy of the system, but also to realize complex dissipators. Generalizing these techniques may allow addressing open questions in non-equilibrium quantum many-body physics.

IX. ACKNOWLEDGMENTS
We thank Achim Rosch, Petar Jurcevic, Christine Maier, and Christian Roos for valuable discussions. We acknowledge financial support from the Swiss National Science Foundation through the National Centre of Competence in Research for Quantum Science and Technology (QSIT) grant 51NF40-160591. FR acknowledges financial support from the Alexander von Humboldt foundation through a Feodor Lynen fellowship and from the Swiss National Science Foundation (Ambizione grant no. PZ00P2 186040). FL acknowledges the financial support of the German Science Foundation under CRC TR 183 (project A01). ZL was financially supported by Gordon and Betty Moore Foundation's EPIC initiative, Grant No. GBMF4545 and from the European Research Council (ERC) synergy UQUAM project, and partially by TR CRC 183 (project A01). ZL also acknowledges the L'Oréal-Unesco national scholarship "For women in science" and the NeTex program of University of Cologne which funded visits of Harvard University where this work was initiated. Here we review the theory of weakly driven and open nearly integrable systems. We consider the setup discussed in the main text, where the dynamics of the system is given by a dominant integrable Hamiltonian H 0 , in the presence of perturbations of unitary and Markovian nature, which weakly break the integrability. The corresponding Liouvillian terms are The steady state density matrix is determined bẏ where L p = L u + L m . Because perturbations due to the Markovian dissipation L m and the next-nearest neighbor interaction L u are only weak, the exact steady state density matrix ρ ∞ can be split as with respect to the eigenstates of H 0 and is parametrized with about 2 N parameters a mn . However, it turns out that description with ρ BD is redundant and can be replaced with a generalized Gibbs ensemble, ρ BD → ρ GGE , if H 0 is integrable [47,48] and with a Gibbs ensemble ρ BD → ρ th , if H 0 is ergodic [87].
While equivalence of ρ BD and ρ GGE is formally expected when calculating expectation values of local observables in the thermodynamic limit and with all conservation laws included, in most cases also a truncated GGE (tGGE) with a few conservation laws qualitatively well captures the expectation values of local observables, if their support is much smaller than the support of included conservation laws. Known exceptions are observables that are orthogonal to all included conservation laws [88].
In the following we will compare the expectation values evaluated with respect to ρ ∞ , ρ BD and ρ tGGE in order to confirm that steady states can be approximately described with non-thermal generalized Gibbs ensembles. Note that since our choice of Lindblad operators breaks magnetization, as well as momentum conservation, the whole Hilbert space is of relevance.
The Lagrange parameters λ i in Eq. (A5) and parameters a mn in Eq. (A4) are determined from the stationarity conditions in the steady state [47][48][49], for ρ tGGE and ρ BD , respectively. For our choice of perturbation the contributions to order and 2 uniquely fix the λ i (or equivalently a mn ) in the steady state. Note that L 0 ρ tGGE = 0 because [H 0 , C i ] = 0, and that unitary perturbation contributes to the decay of conservation laws only in the second order, since tr[C i L u ρ tGGE ] = 0 due to cyclicity of the trace. More details on the derivation of condition (A8) and how to use L −1 0 in practice can be found in Ref. [49]. Here we give only the final result for Ċ i , where finite broadening η has to be used for calculations at finite system sizes. Expressions relevant for ρ BD are obtained by replacements C i → |m n| and ρ tGGE → ρ BD .

Numerical results
We base our analysis on three approaches: (i) calculation of the exact steady state ρ ∞ , Eq. (A1), at finite but small , obtained from diagonalization of the full Liouvillian on small system sizes N = 6 where we exclude or include (NN) unitary integrability breaking H 1 , (ii) exact calculation of ρ BD , Eq. (A4), on N = 6, 8, 10 and (iii) approximate calculation based on a truncated GGE, Eq. (A5), including a finite number of N C = 4 conservation laws C i on N = 10. Note that each C i = i c (i) j is a translationally invariant sum of operators c (i) j with support not larger than i. Due to finite size effects, only C i with support smaller than N/2 can be included in the tGGE. C i are obtained using the so-called boost operator, B = −i j jh j , where H 0 = j h j , from the recursive relation C i+1 = [B, C i ] for i ≥ 2 and C 2 = H 0 . At the isotropic point, J y = J z the magnetization S x = C 1 is conserved as well. Fig. 13 shows H 0 and C 4 as a function of relative dissipator strength γ, Eqs. (10,11), obtained using different approximations described above at largest accessible system sizes. We observe a good agreement between the three approaches, also for other parameters not displayed. Results calculated from ρ BD on N = 6, 8, 10 interpolate between the exact (N = 6) and tGGE (N = 10) results. While ρ BD and ρ ∞ for small = 0.01 agree very well on N = 6, increasing the system size shows a tendency of ρ BD towards the ρ tGGE result. A milder discrepancy of ρ tGGE results is due to omitted conservation laws.
The important conclusion is two-fold: (i) The above analysis gives numerical support for the claim that the stabilized steady state can be approximated with a GGE despite different sources of integrability breaking. (ii) While ρ tGGE is parametrized with N C = 4 parameters, ρ BD at N = 10 with about 10 3 and the full ρ ∞ would require about 10 6 parameters. Description in terms of a truncated GGE is therefore a highly compact parametrization of the steady state, which takes into account only the most relevant information. We find that in the presence of Lindblad driving, the effect of next-nearest interaction, H 1 , is rather weak. While in a closed setup H 1 is crucial as it dictates relaxation towards a thermal state, it is dominated with Lindblad terms in an open setup. Mathematically, this can be explained through Eq. (A10) which shows that the unitary perturbation is constrained to act only between degenerate eigenstates, while the Markovian contribution has no such constraint. Fig. 13 shows that results obtained from the exact steady state ρ ∞ calculated with (NN) or without H 1 are very similar. While H 1 can be easily included into the calculation of the exact steady state ρ ∞ , it brings certain ambiguity into the calculation of ρ BD and ρ tGGE . Namely, on finite system sizes one has to introduce broadening η when calculating L u L −1 0 L u , Eq. (A10). As we showed in [48], broadening itself modifies the effective strength of the perturbation, meaning that different system sizes, requiring different broadening, cannot be directly compared. Since ρ ∞ shows that the effect of H 1 is small, we omit it in the calculation of ρ BD and ρ tGGE .

Appendix B: Microscopic derivation of two-body decay
In the following, we verify that the mechanisms presented in Sec. IV C lead to the desired dissipative couplings in Eq. (11). To this end, we eliminate the excited degrees of freedom by means of the effective operator formalism [76]. This allows us to obtain the effective dynamics of the ground states.
To obtain the effective processes between the ground states, we need to evaluate the expressions for the effective Hamiltonian and Lindblad operators [76], with the relevant terms discussed below.
The jump operators relevant for Eq. (B5) are given by induced spontaneous emission, as described by Eq. (B4). The non-Hermitian terms in Eq. (B5) can then be taken into account by generalizing the detunings from H e to "complex" energies of the form∆ = ∆ − i(Γ/2)/2. Here we assume no motional decoherence and hence,δ = δ. If necessary, processes like phonon decay, L κ = √ κa, can be taken into byδ = δ − iκ/2.
which mediate the effective processes. Using Eq. (B2), we obtain for the effective jump operators for spontaneous emission with the effective decay rate