Dissipative ground-state preparation of a spin chain by a structured environment

We propose a dissipative method to prepare the ground state of the isotropic XY spin Hamiltonian in a transverse field. Our model consists of a spin chain with nearest-neighbour interactions and an additional collective coupling of the spins to a damped harmonic oscillator. The latter provides an effective environment with a Lorentzian spectral density and can be used to drive the chain asymptotically towards its multipartite-entangled ground state at a rate that depends on the degree of non-Markovianity of the evolution. We also present a detailed proposal for the experimental implementation with a chain of trapped ions. The protocol does not require individual addressing, concatenated pulses, or multi-particle jump operators, and is capable of generating the desired target state in small ion chains with very high fidelities.


I. INTRODUCTION
The general goal of quantum information processing is to manipulate the information coded in a particular quantum system, while simultaneously trying to isolate it from undesired perturbations due to the external environment. For instance, in the fields of quantum computation [1] and quantum simulation [2], the processing stage usually relies on a designed unitary evolution which, preserving the coherence of the system, provides a quantum state that contains the outcome of the computation, or the properties of the quantum phase of matter under study. Therefore, the experimental progress in these fields has been typically associated to technological developments that minimize experimental imperfections and maximize environmental isolation [1,2].
However, the coupling of a system to its surrounding environment is not necessarily a disadvantage. In fact, the dissipation of a quantum system, if judiciously exploited, can act as a resource for quantum information processing. This interesting change of paradigm started with the recognition that dissipation can assist the generation of entanglement between distant atoms, either in free space [3], or trapped inside cavities [4,5]. By measuring the presence of a photon spontaneously emitted by the atoms [3,5], or its absence [4], it is possible to project an initially uncorrelated atomic state into a maximally entangled one, as has been recently experimentally demonstrated [6]. These schemes can be improved further by tailoring the atom-cavity interaction [7] or by modifying the measurement techniques [8], and can also be generalized to provide a route towards universal quantum computation [9]. Since the dissipation does not always render the desired result in these approaches, the measurement is required to select only the successful outcomes. Therefore, the combination of dissipation and measurement can be considered as a probabilistic resource for quantum information.
Another promising approach is to assess whether, by designing the system-environment coupling in a particular way, the dissipation could provide the desired quantum state with certainty. This quantum reservoir engineering finds its roots in the theory of laser cooling of atoms [10] and, more related to this work, of trapped ions [11]. Laser-cooling schemes try to design the system-environment coupling such that the dissipation yields a stationary state with reduced kinetic en-ergy. These ideas can be taken a step further in order to nonprobabilistically produce stationary states that display nonclassical aspects, or a certain amount of entanglement [12]. For instance, by mimicking superradiance phenomena [13], certain engineered decay channels yield partial entanglement in the stationary mixed state [14]. However, it would be highly desirable to devise dissipative protocols that provide maximally-entangled pure states asymptotically, the so-called dark states, showing fidelities comparable to those obtained by more standard unitary schemes [1]. This approach has been recently pursued by different groups which have shown that, provided that one can engineer a dissipation that acts quasi-locally on different parts of the system, it is possible to design dissipative protocols to produce a number of paradigmatic dark states [15], or to perform universal dissipative quantum computation [16]. For instance, by engineering a particular dissipation that acts on pairs of adjacent atoms in an optical lattice, the system can be driven dissipatively into a superfluid phase [17]. In addition, as initially proposed for Rydberg atoms [18] and first realized in experiments of trapped ions [19], by concatenating multi-qubit gates with a controlled dissipation of ancillary qubits, a variety of multipartite entangled states (i.e. Bell, Greenberger-Horne-Zeilinger, and Dicke states) have been dissipatively generated with reasonably high fidelities. Let us note that this stroboscopic time-evolution corresponds to a Markovian master equation in the limit of many gates and dissipation steps, where the dissipative jump operators correspond to multiqubit jump operators. From a fundamental point of view, it is of interest to address whether similar dissipative protocols can also work in an (i) analog and (ii) global fashion (i.e. using always-on couplings without individual addressing). Additionally, from a more pragmatic point of view, such global analog schemes would not be limited by the accumulation of the errors in each step of the stroboscopic protocols. However, finding particular schemes to provide multi-qubit jump operators in an analog manner seems to be a tremendous task from both a theoretical and experimental point of view. Therefore, we impose a further constraint, the jump operators should be composed of (iii) single-qubit operators.
In this article, we develop an instance of a global analog dissipative protocol that generates multipartite entangled states corresponding to ground states of a quantum spin chain. The underlying idea is to engineer jump operators that are a particular sum of single-spin operators, in analogy to the models of collective spontaneous emission [13]. We show that, if the dissipation of the spins is mediated by a common harmonic mode, the spin chain sees a structured Markovian environment. By structured environment we refer to a spectral density exhibiting a sharp peak at a certain frequency (so called Breit-Wigner resonance), corresponding to a weakly damped harmonic oscillator. This structure will allow us to design jump operators in such a way that the stationary state of the system corresponds to the ground state of the quantum spin chain. In particular, we consider an isotropic XY spin chain, which describes a critical phase of matter.
We explain in detail how to implement this protocol with trapped ions, relying on tools that have already been achieved experimentally. These tools are the so-called state-dependent forces [20][21][22][23][24], and sympathetic resolved-sideband cooling developed for quantum computation [25,26]. We test the scheme for realistic parameters, and show that the fidelities that can be achieved are comparable to the unitary protocols that produce multipartite entangled states [27]. Our procedure performs well for small chains of trapped ions, which is actually the regime in most of the current experiments.
This article is organized as follows. In Sec. II, we describe the dissipative protocol. We start by introducing the model under consideration in Sec. II A, and then move onto an analytic discussion, supported by numerical results, of how to engineer a structured environment that allows us to prepare multipartite entangled states dissipatively (Sec. II B). In Sec. II C we explain how the method can be used to cool the system to the many-body ground state of the quantum spin chain. In Sec. III we show how crystals of trapped atomic ions are ideally suited for an experimental implementation of our ideas. Section III A describes the trapped-ion setup, a two-species Coulomb crystal. In Secs. III B, III C and III D we introduce the trapped-ion toolbox required to implement the dissipative protocol. Section III E contains a numerical analysis of the trapped-ion dissipative protocol. Finally, we present our conclusions and an outlook in Sec. IV.

II. STEADY-STATE ENTANGLEMENT OF A SPIN CHAIN
A. Spin chain with controllable decoherence Let us start by introducing the model under consideration: an interacting spin- 1 2 chain coupled to a damped bosonic mode (see Fig. 1(a)). This model, which shall be referred to as the damped spin-boson chain (DSBC) in the rest of the manuscript, is a many-body generalization of a system considered recently for the dissipative generation of two-qubit entanglement [28]. The DSBC is described by the following master equation (h = 1) where ρ is the total density matrix, and L DSBC is the dissipative Liouvillian. We consider a finite chain of N spins evolving a J J Figure 1. Damped spin-boson chain: (a) Schematic representation of a spin chain with an isotropic XY Hamiltonian subjected to an additional transverse field whose intensity depends on the position coordinate of a damped harmonic oscillator. This oscillator mediates the dissipation of the spin chain, which effectively sees a structured environment. (b) The effect of the oscillator on the spin chain is described by jump operators L + n and L − n going respectively up and down the spectrum of spin-chain excitations with quasi-momentum q n and energy ε n (the model parameters are defined in the text).
under the isotropic XY Hamiltonian [29] where J > 0 is an antiferromagnetic interaction strength, and we introduced the raising and lowering operators σ The bosonic Hamiltonian consists of a single mode where a † and a are respectively the creation and annihilation operators (the oscillator frequency ∆ a can actually correspond to a detuning ∆ a ≶ 0, as will become clear from the experimental implementation). The spins and the boson are coupled through In the above expression, we have introduced the sitedependent coupling strength g i , and the Pauli matrices Finally, we also consider weak damping of the bosonic mode via a Lindblad-type dissipator where κ > 0 is the damping rate. We note that all these individual ingredients can be realized with the state-of-the-art technology of trapped-ion crystals, as described in Sec. III. The objective of this work is to understand the interplay of these dynamics, such that the Liouvillian L DSBC generates stationary entanglement between the initially uncorrelated spins. So far, we can already appreciate two of the properties of the protocol outlined in the introduction: (i) It is analog, since the couplings (2)-(5) are considered to be switched on during the whole protocol. (ii) It is global, since the couplings (2)-(5) address all the spins of the chain. Let us note that the spin-chain Hamiltonian (2) alone could already generate entangled states unitarily at given instants of time. Nonetheless, these correlations have a transient nature, whereas the interest of this work is the onset of stationary correlations. Therefore, we explore the interplay of this Hamiltonian part with the irreversible dynamics introduced by the damped boson through the spin-boson coupling (4).
In Fig. 1(a), we represent schematically the DSBC. The spin chain is considered to be our system S , whereas the bosonic mode, together with the Markovian reservoir where it dissipates, provide an effective structured environment E for the spin chain. At first sight, the dephasing-like spin-boson coupling (4) seems to introduce a source of decoherence in the spin chain, thus hindering rather than assisting the generation of stationary entanglement. In the following section, we will show that this naïve intuition is not always valid, and explain the subtle mechanism that allows for generation of stationary entanglement in a particular regime of the system parameters.

B. Cooling or heating to an entangled steady state
We consider an initial state of the system is an arbitrary state of the harmonic mode, and the uncorrelated spin state |ψ 0 = |↑ 1 ↓ 2 ↓ 3 · · · ↓ N has a single spin excitation at the left edge of the chain. In this section, we discuss how the system Liouvillian L DSBC acts on this state, allowing for the onset of stationary entanglement.

Spin-chain spectrum
As a consequence of the XY Hamiltonian (2), the initially localized spin excitation is exchanged between different neighbors. In fact, the spin chain corresponds exactly to a tight-binding model where the excitation hops along the sites of the underlying chain according to In this expression, H s,1 results from the projection of the XY Hamiltonian onto the single-excitation subspace H s → P s H s P s , where P s is the projector onto an N-dimensional subspace spanned by the states |i , such that i labels the position of the spin excitation in the chain. This tight-binding model gives rise to the following band structure (see Fig. 1 where each energy is associated to a "spin-wave excitation" The goal of this section is to find the conditions to generate dissipatively one of these multipartite entangled states, Tr a {e L DSBC t |ψ 0 ψ 0 | ⊗ ρ b (0)} → |q n q n |.

Spin-wave ladder
We consider a spatial modulation of the spin-boson coupling strength (4) in the form g i = g cos(q g i), where q g = π/(N + 1). In the single-excitation subspace, the spin-boson coupling becomes where the ladder operators L + n = |q n q n+1 | = (L − n ) † are responsible for ascending or descending along the ladder of spin-wave excitations (see Fig. 1(b)). Accordingly, the spinboson coupling (9) connects different spin-wave excitations, while simultaneously creating or annihilating bosonic quanta. Moreover, since the harmonic mode dissipates irreversibly into a Markovian environment (5), the ladder operators L ± n introduce irreversible dynamics in the spin chain. In terms of these ladder operators, the single-excitation spin-chain Hamiltonian (6) reads where for simplicity we reset the zero of energy at the bottom of the band. We show below how to generate mesoscopic entangled states by controlling the relative strengths of the dissipative processes that take the system up and down the spinwave ladder.

Irreversible dynamics in the spin chain
Since we are interested in the stationary properties of the DSBC, we consider long times t κ −1 , such that the bosonic mode has enough time to relax under the action of the damping. Besides, for |g| κ, this relaxation is much faster than the process of energy exchange between the boson and the spins. In this limit, the boson thermalizes individually, and the spin chain evolves on a much slower timescale under such a bosonic background.
To obtain the effect of the boson background on the slower spin dynamics [30], we must "integrate out" the bosonic degrees of freedom from the spin-boson coupling (4). We start with a state of the form ρ(t) = ρ ss b ⊗ ρ s (t), where ρ ss b is the vacuum of the harmonic mode, and fulfills D b ρ ss b = 0. We then expand the state at a later time t + δt, with κδt 1, in powers of the coupling constant g between the spin system and the oscillator, keeping up to second order in g. Tracing over the mode, we obtain: where the "hat" indicates that we work in interaction picture with respect to H s,1 + H b , and where we have introduced Using the explicit form of D b and L sb , after some algebra and one integral in t we find: Here, we introduced the collective jump operators with ξ ± n = κ/[κ + i(∆ a ± ∆ n )], being ∆ n = ε n − ε n+1 > 0 the energy difference between two neighboring spin-wave excitations in the spin-wave ladder (see Fig. 1(b)).
The integrand in Eq. (12) contains some terms that oscillate rapidly in time. We now perform the remaining integral, assuming that the time interval δt is short compared to the time scale given by κ/g 2 , but long enough such that δt|∆ n −∆ n | 1 ∀∆ n = ∆ n . This condition restricts the values of g for which this treatment is valid. To perform the integral, we group the frequencies ∆ n in such a way that within each group the frequencies are equal, and the difference between the frequencies in different groups is finite (we note that this is not possible in the limit of infinite sizes, where the energies become a continuum). It is worth noticing that the presence of degenerate frequencies is typical in this model and makes the grouping necessary. If one keeps only the dominant terms, which are the terms in the integrand that are constant in time, back in Schrödinger picture one finds the following master equation governing the coarse-grained evolution over the time scales of interest: with and where In the expressions above, the sum over ∆ runs over the different transition frequencies in the system, { , } denotes an anticommutator, and the Lindblad operators are defined as: where the sum is over all the values of n such that ∆ n = ∆. The action of the different Lindblad operators is weighted by the spectral density of the effective environment, According to these expressions, the spins are subjected to a Lorentzian reservoir centered at the boson detuning ∆ a with a width given by the boson damping rate κ. Therefore, the dissipation on the spins is not equal at all frequencies, but stronger at frequencies matching that of the bosonic mode (i.e. they do not see a totally flat environment, but a structured one [31]). As announced previously, the ladder operators L ± n are responsible for introducing the irreversible dynamics in the spin chain. In particular, they determine the collective jump operators (17), where the adjective collective emphasizes that they act over all the spins in the chain. However, we remark that these jump operators are a sum of single-spin operators, as opposed to the multi-spin nature of some other engineered dissipation protocols considered recently [15,16,18,19]. With this discussion, we show the third property of the protocol outlined in the introduction: jump operators are composed of sums of (iii) single-qubit operators. This draws an analogy to the models of collective spontaneous emission [13], but we will show that the special dissipation mediated by the boson mode allows us to use them to prepare dissipatively the ground state of the spin chain.

Stationary entanglement
We are searching for the steady state ρ ss s of Eq. (14), L s (ρ ss s ) = 0, with three properties: being pure, displaying multipartite entanglement, and being unique. A pure steady state ρ ss s = |Ψ ss Ψ ss | necessarily belongs to a Hamiltonian eigenspace [15], which in our case is given by the spinwave excitations |Ψ ss ∈ {|q n , n = 1 · · · N}. Additionally, the steady state should also belong to the kernel of the dissipator, D s (|Ψ ss Ψ ss |) = 0. By inspection of the jump operators, we identify two regimes where this can happen, which depend on the boson detuning with respect to the spin-wave energy difference: For the regime of positive detunings, we can then approximate the evolution using that J a (−∆ n ) is negligibly small for all n, while all the J a (∆ n ) are non-negligible. In that case the nonunitary part of the evolution is taking the system to the lowest state in the manifold of spin waves (the Hamiltonian part does not induce transitions between spin waves). Thus, the unique pure steady state corresponds to the lowest spin wave |Ψ ss = |q N . Conversely, for Figure 2. Tailoring the dissipative jump operators: (a,c) In the regime of positive oscillator detuning ∆ a ≈ ∆ n , the jump operators that descend the ladder govern the dissipative dynamics (i.e. r n 1). This can be understood in terms of the effective Lorentzian spectral density J a of the environment seen by the spin chain, which overlaps minimally with the processes that climb up the ladder. (b,d) In the regime of negative oscillator detuning ∆ a ≈ −∆ n , the jump operators that climb up the ladder dominate the dissipative dynamics (i.e. r n 1), since in this case the spectral density of the environment seen by the spin chain overlaps minimally with the processes that descend the ladder. negative detunings, the steady state corresponds to the highest spin-wave |Ψ ss = |q 1 . Following [15], it can be shown that in these limits the particular form of the Lindblad operators does not allow for mixed steady states. We thus get the desired result: a unique pure state displaying stationary multipartite entanglement.
In Fig. 2, we represent schematically the two regimes of interest. When conditions (19) are fulfilled, for κ ∆ a ≈ ∆ n , the dissipative processes that go down the spin-wave ladder dominate since |ξ + n | |ξ − n |. Effectively, the dissipation induced by the damped bosonic mode "cools" the spin state to the lowest-energy spin wave, Fig. 2(a)). The bath spectral density peaks at the frequencies corresponding to the dissipative processes descending the spin-wave ladder (see Fig. 2(c)) and therefore, the bath absorbs very efficiently the energy dissipated by the spin chain during its relaxation to the spin wave |Ψ ss = |q N . Conversely, when the system parameters fulfill conditions (20), for negative detunings κ −∆ a ≈ ∆ n , it is the dissipative processes that climb up the ladder which dominate, |ξ − n | |ξ + n |. In this limit, the bosonic mode drives the spin state to the highest-energy spin wave, Fig. 2(b)). The bath spectral density is maximal for the ascending dissipative processes (see Fig. 2(d)), and the bath provides the required energy to climb up the spin-wave ladder. We may thus conclude that in both regimes, the structured reservoir singles out a unique spin wave as the steady state of the chain, assisting the generation of stationary multipartite entanglement.
In the light of these results, we can also understand the regime κ J. In this limit, the reservoir spectral density (18) becomes essentially flat, and both up/down processes contribute equally ξ + n ≈ ξ − n ≈ 1. It is straightforward to see that if J a (ω) is constant, the dissipator in Eq. (16) satisfies D s (I) = 0, and the totally mixed state ρ ss s ∝ I becomes a steady state. This regime corresponds to the naïve argument of section II A, following which one expects that the dephasing-like term (4) can only decohere the spin chain. It is now clear that there are other regimes, (19)- (20), where we can profit from a structured environment to assist the generation of stationary entanglement.
We now comment on the time required to reach the aforementioned steady states. In the regime where the effective Liouvillian (14) was derived, it is given by For a particular experimental setup, where the spin couplings cannot reach arbitrarily large values, this preparation time may turn out to be too long for practical purposes. We emphasize, however, that the same dissipative preparation of entangled states can be obtained in a non-Markovian regime where |g| ∼ κ J, or κ |g| J. Since simple master equations cannot be derived analytically for this regime, we shall explore it numerically.

Dissipative preparation of W-like states
Let us note that although the spin-wave excitations (8) are genuinely multipartite entangled, the weights for the excitation of each of the spins in the system are different. A small modification of the scheme also allows for the dissipative generation of N-partite W-like states, where the W-state is defined as Let us leave for a moment the constraint to consider purely global couplings, and assume that it is possible to add external transverse fields that act locally at the two edges of the chain In the subspace with one spin excitation, the ground state of this Hamiltonian is of the form which is locally equivalent to the above N-partite W-state (we note that the actual W-state can be obtained if the spin-spin coupling is ferromagnetic instead of antiferromagnetic). As opposed to the ground state of the spin system considered so far, the W-state has the property that it is fully symmetric under particle interchanges. In the state |q N , the excitation is equally distributed over all the sites, and each of the particles is equally entangled with the rest. This state can be prepared with the same method described before, at the price of introducing individual addressability in the trapped-ion proposal of Sec. III.

Numerical analysis of the protocol
In this section, we analyze numerically the validity of our previous analytical derivations by integrating directly the master equation of the damped spin-boson chain (1). We obtain the fidelity with the desired stationary entangled state where |ψ 0 = |↑ 1 ↓ 2 ↓ 3 · · · ↓ N , we consider that the boson mode is initially in the vacuum state ρ b (0) = |0 0|, and |Ψ target is the particular entangled state that the protocol targets. We test the robustness of this fidelity against a nonoptimal choice of the system parameters (19)-(20), a limited protocol time t f , and an increasing particle number N.
In Figure 3, we consider the time evolution of the fidelity for the preparation of the lowest-energy spin-wave excitation of an N = 3 spin chain For these calculations, we set the oscillator detuning ∆ a to the value ∆ a = √ 2J which is resonant with the spin-wave transitions, and vary the ratio between the spin-boson coupling strength g and the damping rate κ. The plot indicates that, for similar values of the final fidelity, the convergence is much faster when g and κ are of the same order. This implies that for practical purposes it is most convenient to work in the deeply non-Markovian regime. The optimal ratio g/κ, however, depends in general on the number of spins and the value of ∆ a .
It is important to note that the resonance condition can be relaxed to some extent. In Figure 4 we analyze the fidelity as a function of the detuning ∆ a and of the coupling g, setting κ = g. From our previous analysis, we expect the fidelity to be maximal when ∆ a = ∆ a , and when g, κ approach zero to fulfill the requirement g, κ J. This behavior is totally captured by the numerical results in Fig. 4(a), displaying the asymptotic fidelity. The maximum fidelity, highlighted with a white dot, is indeed found in this limit and is equal to 1 within the numerical errors. We remark that, even when the above conditions are only partially fulfilled, the dynamics can provide the desired state with fidelities well above 90%. Therefore, we can claim that our dissipative protocol is considerably robust against parameter imperfections.
Nevertheless, the time required to approach the desired target state might be prohibitively long. Therefore, we also characterize the fidelity for a fixed protocol duration. In Fig. 4(b), we show the value of the fidelity at a finite fixed time t f = 10 3 /J. As can be seen in this figure, the optimal choice for the oscillator detuning is still ∆ a = √ 2J, while the optimal value for g = κ has moved away from zero. In fact, the best choice now results from the interplay between the condition g, κ ∆ a , J, and the required convergence speed. We emphasize that it is still possible to find parameters such that the achieved fidelities remain very high, F t f = 0.9999.
We now address the dissipative generation of W-like states by adding boundary transverse fields (23). For the N = 3 spin chain under consideration, the target state is In Fig. 4(c)-(d), we show the asymptotic and finite-time fidelities for the dissipative generation of such a tripartite entangled state with homogeneously-distributed excitations. We observe high fidelities comparable to the inhomogeneous spin-wave state (8), and moreover, an analogous robustness with respect to a non-optimal choice of the system parameters. Let us note that, for N = 3, there are two different transitions in the spectrum of these W-type spin-waves: one is resonant at a frequency of∆ 2 = J, and the remaining one occurs at∆ 1 = 2J. Therefore, it is impossible to set an oscillator detuning at resonance with both transitions simultaneously, and thus the conditions (19) are not completely fulfilled. This explains the fact that the W-type fidelities achieved for finite times are lower than those corresponding to the inhomogeneous spin waves in Fig. 4(b), where the two transitions have the same frequency In this section, we have analyzed numerically the validity of the scheme for the generation of the simplest case of multipartite entanglement, namely tripartite entangled states. A question that should be carefully addressed is whether the same scheme can be used to generate N-partite entangled states, and how large can the attained fidelities be as N is increased. This is the topic of the following section.  (19) is met for ∆ a = J, and g J). (d) The same as in (c) but setting a finite time t f = 10 3 /J. For these numerical calculations, the oscillator was truncated to three levels, the trace is preserved to 10 −6 and the maximum fidelity is equal to 1 within this error. In each plot, the highest fidelity is indicated by a white dot.

Mesoscopic spin chains
We start by assessing the fulfillment of the necessary conditions (19)- (20), which rely on the energetic difference between the dissipative processes that climb up/down the spinwave ladder, as the number of spins N is increased. For very large spin chains N → ∞, the energy differences between neighboring spin waves |∆ n | 2πJ/N → 0. This implies that the energetic argument selecting only processes going up or down the ladder can no longer hold. Indeed, lim N→∞ r n = lim N→∞ |ξ + n |/|ξ − n | = 1, so that the protocol ceases to be operative in the thermodynamic limit of infinitely long chains. However, for mesoscopic spin chains, the ratio r n can be controlled to an acceptable degree. In Fig. 5(a), we show that for positive detunings and N ∼ 10, r n ≈ 0.1 for the most of the spin-wave excitations, whereas slightly higher values are attained for the extremal spin excitations (where |∆ n | ∝ J/N 2 ). For ratios r n on this order, we expect that the dissipative preparation of entangled states still works with acceptable fidelities.
To be more specific, we note that the presence of different transition frequencies ∆ n is generic for N ≥ 4, which represents an obstacle for the efficiency of the procedure. As the number of spins is increased, the convergence also becomes slower because of the larger number of steps down/up the ladder towards the target state, and the lower transition frequencies which require a decrease in g and κ. In Fig. 5(b), we show numerical results for the dependence of the protocol error ε t f achieved for a finite time t f = 10 3 /J as a function of the number of sites and optimizing the values of g, κ and ∆ a .
We observe that the errors obtained for the dissipative state preparation of the inhomogeneous spin wave |q N are below 10% for chains of up to ten sites. In comparison, the creation of the W-like states is worse, and such high fidelities can only be achieved for short chains of up to five sites.
C. Preparation of the ground state of the isotropic XY chain in a transverse field So far, our analysis has been restricted to an N-dimensional subspace of the spin chain, since the initial spin state |ψ 0 = |↑ 1 ↓ 2 ↓ 3 · · · ↓ N contains a single excitation, and the number of spin excitations is preserved by the complete DSBC Liouvillian (1). In the following we explore the full 2 N -dimensional Hilbert space of the spin chain.

Jordan-Wigner jump operators
In order to treat the full Hilbert space of the spins, we fermionize the spin chain via the so-called Jordan-Wigner transformation [32], namely where c † i , c i are fermionic creation-annihilation operators. The spin-chain (2) and spin-boson (4) Hamiltonians can be Error ε t f = 1 − F t f for the dissipative generation of the target state in a finite time t f = 10 3 /J as a function of the number of spins N in the chain. The blue circles represent the errors corresponding to the inhomogeneous spin waves |q N ∝ sin(q N )|↑↓↓ · · · ↓ + sin(2q N )|↓↑↓ · · · ↓ + · · · + sin(Nq N )|↓↓↓ · · · ↑ . The red diamonds represent the error for the W-like states |q N ∝ |↑↓↓ · · · ↓ − |↓↑↓ · · · ↓ + · · · + (−1) N+1 |↓↓↓ · · · ↑ . For the numerical calculation, the Hilbert space of the oscillator was truncated to four levels. expressed in terms of the Jordan-Wigner fermions as follows where we recall that the spin-boson couplings are g i = g cos(q g i) with q g = π/(N + 1). The next step is to express these Hamiltonian terms in the spin-wave basis introduced in Eq. (8), which leads to the following expressions where the fermionic operators in momentum space are c q n = N ∑ i sin(q n i)c i . In combination with the dissipative part given by Eq. (5), the spin-boson system is analogous to a damped single-mode Holstein model, a dissipative version of the familiar Holstein model describing electron-phonon interactions [33].
We can obtain the same formal expressions as in the single-excitation problem by rewriting the ladder operators in second-quantized form Accordingly, in the regime where the boson degrees of freedom can be integrated out, we obtain a purely fermionic master equation which coincides with Eqs. (14)- (16), but with the collective jump operators (17) now expressed in terms of the Jordan-Wigner ladder operators As a result, the dissipative dynamics restricted to the singleexcitation sector can be generalized directly to the full Hilbert space with arbitrary numbers of spin excitations.

Effective ground-state cooling
Let us now consider an initial state with an arbitrary number of spin excitations n s ≤ N distributed along the chain We would like to determine the steady-state of the spin-boson system if the conditions (19) are fulfilled. In this limit, the Lindblad operators are only of the form Such a jump operator has the effect of lowering the energy of the fermionic quasiparticles. Due to the Pauli exclusion principle, the initial n s excitations cannot all occupy the lowest energy level, with quasimomentum q N . Instead, the stationary state must be of the following form This is precisely the ground state of the original isotropic XY model (2), if supplemented with a homogeneous transverse field H s → H(J, µ) = H s − (µ/2) ∑ i σ z i . In this Hamiltonian, the transverse field plays the role of an effective chemical potential µ = 1 2 (ε N−(n s −1) + ε N−n s ), which is determined by the initial number of spin excitations. Due to the to the dissipative process, the excitations are distributed in the lowest available single-particle states (see Fig. 6(a)). This means, for each number of spins up in the initial state prepared, there is a value of the transverse field µ such that the asymptotic state corresponds to the ground state of H(J, µ). Conversely, a given choice of J, µ determines the number of spins that should be up in the initial state so that the dissipative dynamics take the system into the desired ground state. Figure 6. Cooling to the ground state of the spin chain: (a) Schematic representation of the dissipative processes that cool each of the n s spin excitations to the lowest non-occupied state. The stationary state corresponds to the ground state of the isotropic XY Hamiltonian, with the particular filling determined by an additional transverse field acting as a chemical potential. (b) The same as in (a) but selecting the processes that climb up the ladder. The stationary state would correspond to the zero-temperature state of the same XY model in a transverse field, but with reversed coupling strengths.
We could also explore the steady state when the conditions (20) are fulfilled. In this limit, the only Lindblad operators are of the form: and pump all the excitations to the highest-energy available single-particle states corresponding to the ground-state of the XY Hamiltonian H(−J, −µ) with a ferromagnetic spin-spin coupling, and an inverted chemical potential (see Fig. 6(b)).

Numerical results
Let us now confirm the validity of the above argument by numerical analysis for small spin chains with different numbers of initial excitations. In Fig. 7, we display the error in producing the ground-state (35) taking a fixed protocol time of t f = 10 3 /J. The highest errors occur when the number of excitations is 1 or N − 1, since the transition frequencies are lowest at these points (this is the worst case, with the errors plotted in Fig. 7 by means of red squares). With yellow squares, we represent the lowest attained errors. Note that the average errors of the protocol (blue circles, averaged over all the possible excitation numbers) are as low as 10 −2 for chains up to N = 6 spins, which would thus allow us to prepare dissipatively the ground state of the spin model with fidelities of 99%.
Let us close this section by noting that, in the general case with an arbitrary number of spin excitations, a different choice for the site-dependence of the coupling coefficients g i may allow for jump operators that go down/up the spin ladder in larger steps, and therefore favor a faster convergence to the asymptotic state (this appears to be specially well-suited for Scaling of the ground-state-cooling protocol: Error ε t f in producing the desired ground-state |G s after a fixed time t f = 10 3 /J and for a number of spins ranging from N = 2 to N = 6. For each number of sites, we plot in blue circles the average over the different possible numbers of excitations (n s ∈ [1, N − 1]), and with yellow and red squares the best and worst cases, respectively. For the numerical calculation, the Hilbert space of the oscillator was truncated to three levels.
chains around the half-filling condition). However, for up to N = 6 spins, we have not found any sizable advantage, and the highest fidelities were always achieved for transitions between neighboring spin waves n → n ± 1.

A. Trapped-ion crystal as a spin-boson chain
Let us describe a specific platform where the DSBC model (1) can be realized. We consider a system of N + N ions confined in a linear radio-frequency trap [34], and arranged forming a one-dimensional chain (see Fig. 8(a)). The vibrational degrees of freedom around the equilibrium positions can be expressed in terms of the so-called normal modes [35], namely where α ∈ {x, y, z} refers to the trap axis, n ∈ {1, . . . , N + N } labels the different modes with frequencies ω nα , and b nα , b † nα , are bosonic operators that annihilate/create phonons in a particular mode. In a linear chain of ions with equal mass, the phonon branches have the structure shown in Fig. 8(b), such that the different modes spread around the trap frequencies lying in the range ω x /2π, ≈1-10MHz, ω y /2π, ≈1-10MHz, and ω z /2π ≈0.1-1MHz. A property of the longitudinal branch that will be used in this work is the existence of an energy gap between the lowest-energy mode (i.e. the center-of-mass mode) and the following one.
As shown in Fig. 8(a), N ions of the chain correspond to a particular atomic species with hyperfine structure, whereas the remaining N ions, that will be used for cooling, do not necessarily have hyperfine splitting (we note that though we . Trapped-ion implementation: (a) Mixed Coulomb crystal in a linear radio-frequency trap. We consider two different types of ions with masses m,m, either of a different species m =m or a different isotope m ≈m. The crucial property is that each type has a very different transition frequency ω 0 =ω 0 . Colored arrows represent the combination of beams in a different configuration (i.e. R= Raman, tw=traveling wave, sw=standing wave) leading to the desired laser-ion interaction. Let us note that the standing-wave state-dependent forces can also be replaced by traveling waves without compromising severely the fidelities of the dissipative scheme. (b) Scheme of the different phonon branches for a two-isotope Coulomb crystal. The vibrational frequencies ω nα span around the trap frequency of each axis ω α . We note that the energy gap between the longitudinal center-of-mass mode and its neighboring mode is equal to ( √ 3 − 1)ω z for the case of ions with equal masses. (c) Atomic Λ-scheme for the N ions with hyperfine structure. A couple of laser beams with Rabi frequencies Ω L α,1 , Ω L α,2 connect the two hyperfine states {|↑ , |↓ } to an auxiliary excited state |aux . When ω L α ≈ ω 0 , the lasers drive a two-photon Raman transition, and a pair of Raman beams of this type lead to the state-dependent forces in Eq. (43). When ω L α ≈ ω nz ω 0 , the lasers lead to a differential ac-Stark shift, such that the contribution from the crossed beams leads to the state-dependent force in Eq. (49). (d) Atomic scheme for the (N + 1)-th ion. A driving of the transition with Rabi frequencyΩ L α is red-detuned from the atomic transitionω L α =ω 0 +∆ L α , such that ∆ L α < 0, leading to an effective damping of the modes. focus on hyperfine qubits, our proposal could be generalized to optical or Zeeman qubits). For simplicity, in the following we describe the case N = 1. For the first N ions, we select two hyperfine levels of the ground-state manifold, which are referred to as spin states {|↑ i , |↓ i }, and have a transition frequency ω 0 in the microwave range with a negligible linewidth Γ ≈ 0 (Fig. 8(c)). For the remaining (N + 1)-th ion, we select two internal states {|g , |e } from a certain transition with frequencyω 0 , and linewidthΓ/2π ≈10-100 kHz (Fig. 8(d)). This is the so-called resolved sideband limitΓ ω z , which will allow for an efficient laser cooling close to the vibrational ground state for the longitudinal center-of-mass mode [25]. Let us note that these small linewidths can be obtained by a variety of methods depending on the particular ion species (e.g. for weakly-allowed dipole transitions, for quadrupoleallowed transitions with an additional laser that admixes the excited state with that of a dipole-allowed transition, or by a Raman configuration) [36].
The complete degrees of freedom are thus described by the following master equatioṅ where we have introduced the Hamiltonians such that τ z i = |↑ i ↑ i | − |↓ i ↓ i |, andτ z N+1 = |e e| − |g g| (we use τ instead of σ to emphasize the pseudo-spin character of these internal states). Additionally, we have to take into account the dissipation from the dipole-allowed transition considering recoil effects [37]. For trapped ions, this can be described by the following super-operatorD(ρ) = D 0 (ρ) +D 1 (ρ), wherẽ whereτ + N+1 = |e g| = (τ − N+1 ) † . This is the usual dissipator of a two-level atom in an electromagnetic environment.
The spontaneous emission also has effects on the vibrational degrees of freedom due to the photon recoil. To leading order in the Lamb-Dicke parameter η nα = (ω 0 /c)/ √ 2mω nα for each normal mode, the jump operators create and annihilate phonons as given by [38]: where we have introducedΓ nα ∝Γ αη 2 nα (M α N+1,n ) 2 , with M α N+1,n the coefficients describing the displacement of the N + 1 ion in mode n along direction α.
Once the vibrational and atomic degrees of freedom have been introduced, let us summarize the role that each of them will play in the realization of the damped spin-boson chain (1): (i) The hyperfine levels of the N ions will simulate the spins in the DSBC model. (ii) The 2(N + 1) vibrational modes along the x and y axes will be used to mediate the spinspin interactions leading to the isotropic XY model in Eq. (2). (iii) The longitudinal center-of-mass mode will play the role of the harmonic oscillator in the DSBC model. (iv) The atomic levels of the (N +1)-th ion will be used to laser-cool the longitudinal center-of-mass mode. We note that the position of this auxiliary ion is irrelevant for the purpose of cooling, but it modifies the spin couplings in the chain. Moreover, it is also possible to use more than one auxiliary ion as an additional knob to control the damping rate. In the following sections, we will detail the different elements required for the ion-trap implementation of the DSBC. We emphasize that all of them can be realized with state-of-the-art technology.

B. Engineering the isotropic XY spin model
We start by addressing how to implement the isotropic XY model (2) in the ion chain by means of a variation of the Mølmer-Sørensen gate [39]. The idea is to use the so-called spin-dependent forces, a tool that has already been implemented in several laboratories for different purposes [20][21][22][23][24]. By combining a pair of Raman laser beams (see Fig. 8(c)), each of which is tuned to the red/blue first vibrational sideband of the hyperfine transition [40], it is possible to obtain the following state-dependent forcê where the "hat" indicates the interaction picture with respect to H 0 = H σ +H σ + H p . In the above expression, we have assumed that the effective wavevector of the Raman beams k L α = k L α,1 − k L α,2 points along the α = {x, y} axis of the trap (see Fig. 8(a)). Additionally, we assume that the strength of both Raman beams is equal, such that they share a common Rabi frequency |Ω L α |, and have opposite detunings, such that we can define a common δ nα = ω L α − (ω 0 − ω nα ).
We have introduced the laser Lamb-Dicke parameter η nα = e α · k L α / √ 2mω nα 1. The constraints on these parameters are |Ω L α | ω α to neglect other terms in addition to such a state-dependent force. Finally, we have defined the sum and difference of the Raman beam phases and the Pauli spin operator τ In order to obtain the desired isotropic XY model, we apply two orthogonal state-dependent forces with the following phases and directions Altogether, the Hamiltonian of the ion chain in the interaction picture becomeŝ Performing a Magnus expansion [41], the evolution turns out to be given by the unitary operatorÛ(t) = exp{Ω eff (t)}, where Ω eff (t) displays the following form to second order We want to derive an effective Hamiltonian from this expression Ω eff (t) ≈ −itH ti s (where the superscript "ti" stands for the trapped-ions realization). The first-order contribution leads to a couple of orthogonal state-dependent displacements, which can only be neglected in the limit |F α in | δ nα [42]. In addition, from the non-commutativity of the σ x and σ y forces, we get a residual spin-phonon coupling in the second-order term of the Magnus expansion. In order to neglect it, we have to impose a further constraint, namely |F x in (F y im ) * | |δ nx − δ my |, which can be fulfilled if the trap frequencies along the x, y axes are sufficiently different ω x = ω y . Under these constraints, an interacting quantum spin chain is obtained (48) Therefore, by adjusting the strengths of the Raman beams, and their detunings, it is possible to find a regime where J x i j = J y i j , and the above spin Hamiltonian corresponds to the desired isotropic XY model since τ . Let us remark that the resulting XY model has the peculiarity of displaying long-range couplings. In fact, when ω x , ω y ω z , and the Raman lasers are far-detuned from the whole vibrational branch, it can be shown that these couplings decay with a dipolar law. The trapped-ion Hamiltonian (48) becomes a realization of the XY interaction (2) in the DSBC model (1) after the identifications τ ± i ↔ σ ± i , and J x i j = J y i j ↔ J. For typical nearest-neighbor distances of z 0 ≈1-10µm, the spin-spin couplings attain strengths in the J α ii+1 /2π ≈1-10 kHz. Let us note that modifications of the scheme of state-dependent forces have been proposed to yield other quantum spin models [42,43], some of which have also been realized experimentally [23].

C. Controlling the spin-boson coupling
Let us now address how to implement the spin-boson coupling (4) in the ion chain. The idea is again to use a Λbeam configuration (see Fig. 8(c)), but in a different regime. Rather than tuning the two-photon frequencies ω L α to the vibrational sidebands of the hyperfine transition, we impose that ω L α ≈ ω nz ω 0 . Moreover, we consider that the laser beams form a standing wave along the trap z-axis (see Fig. 8(a)). Under these constraints, the laser-ion interaction leads to a crossed-beam differential ac-Stark shift, which can be interpreted as another state-dependent force in the τ z basiŝ where the "hat" refers to the interaction picture with respect to H 0 = H τ +H τ + H p . Here, the parameters are defined in analogy to those in Eq. (43), with three important differences: (i) Ω L z is not the two-photon Rabi frequency of the hyperfine transition, but rather a differential ac-Stark shift coming from processes where a photon is exchanged between the pair of laser beams in the Λ configuration. (ii) The detunings are changed to δ nz = ω L z − ω nz . (iii) Since the ion chain lies along the z axis, there is an additional site-dependent phase when shining the lasers such that k L z · r 0 i = 0. When adjusting the difference of the laser phases to be ϕ = π/2, then we obtain F z in = 1 2 η nz Ω L z M z in cos(k L z · r 0 i ). Let us highlight that the standing-wave nature of this spin-dependent force is not essential for the dissipative protocol. However, it makes the connection with the DSBC studied in Sec. II more transparent. We will show in Sec. III E that essentially the same fidelities can be achieved for a traveling-wave configuration, which has also been experimentally demonstrated [40].
We now exploit the particular properties of the longitudinal phonon branch (see Fig. 8(b)). More precisely, we use the presence of a gap between the axial center-of-mass mode and the rest of the modes along the same direction. For the case of equal ions, |ω nz − ω 1z | ≥ ( √ 3 − 1)ω z , while for different ion species or isotopes, the exact value of the gap will depend on the mass ratio and the position of the cooling ions. We assume that ω L z is red-detuned with respect to the axial centerof-mass mode ω L z = ω 1z − δ 1z , such that δ 1z ω 1z . Since the next vibrational modes are separated by a large energy gap, the laser coupling to the remaining phonon branch is highly offresonant and can be thus neglected. Accordingly, this statedependent force giveŝ We can finally move to a different picture where the above Hamiltonian becomes time-independent H L z = H ti b + H ti sb , where Therefore, it is clear that the center-of-mass mode plays the role of the harmonic oscillator {b 1z , b † 1z } ↔ {a, a † } in the DSBC, with the laser detuning corresponding to the oscillator detuning δ 1z ↔ ∆ a in the DSBC model (3). In addition, the strength of the state-dependent force determines the spinboson coupling F z in ↔ g i in the DSBC model (4). Let us now comment on realistic values for the trappedion parameters. Since the laser detuning should only fulfill δ 1z /2π ω 1z /2π = ω z /2π ≈1-10 MHz, it will be easy to reach the required condition of δ 1z /2π ∼ J α ii+1 /2π ≈1-10 kHz. On the other hand, we know from the previous section that g J is necessary to have an accurate pumping to the desired entangled state. Nevertheless, g should not be too small that the total preparation time becomes prohibitively long. A suitable choice could be g ≈ 0.1-0.5 kHz. Since F z i1 = η 1z Ω L z cos(k L z · r 0 i )/ √ N + 1, it will suffice to set the Rabi frequency in the Ω L z ≈1-5 √ N kHz. Finally, we should adjust the laser wavevector such that k L z · r 0 i ≈ πi/(N + 1).

D. Effective damping of the bosonic mode
Once the implementation of the XY model (2) and the spinboson coupling (4) has been described, let us turn into the last required ingredient: an effective damping of the bosonic mode (5). As discussed previously, the idea is to exploit the atomic levels of the (N + 1)-th ion to laser-cool the longitudinal center-of-mass mode. Since the resonance frequencies are very different, ω 0 =ω 0 , the cooling lasers do not affect the spin dynamics of the other N ions. However, since the vibrational modes are collective, it is possible to sympathetically cool the vibrations of the crystal by only acting on the (N + 1)-th ion. This sympathetic laser cooling [44] has been already realized experimentally in small crystals for quantum computation [25,26].
Our starting point is the master equation (39) with the dissipation superoperators (41)- (42). To control the effective damping, we introduce a laser beam red detuned from the atomic transitionωL z ≈ω 0 +∆L z , such that∆L z < 0 (see Fig. 8(d)). Let us note that the particular laser-beam configuration will depend on the particular ion species, and the scheme to attain the resolved-sideband limit [36]. We consider that the laser is in a traveling-wave configuration, such that its wavevector is aligned parallel to the trap axisk L z e z (see Fig. 8(a)). After expanding in series of the Lamb-Dicke parameterη nz = e z ·kL z / √ 2mω nz 1, the laser-ion interaction leads to two different terms, the so-called carrier term and the red and blue sideband terms where we introducedF z N+1,n = i 1 2η nzΩL z M z N+1,n e ikL z ·r 0 N+1 , and the corresponding Rabi frequencyΩL z . By controlling appropriately these sideband terms, and their interplay with the atomic spontaneous emission, it is possible to tailor the damping of the longitudinal modes.
Let us rearrange the full Liouvillian (39) as a sum of two terms L = L 0 + L 1 , where which is justified for small Lamb-Dicke parameters. To obtain the effective damping of the longitudinal center-of-mass mode, we use a similar formalism as in [30]. In this case, the fastest timescale in the problem is given by the decay ratẽ Γ of the atomic states of the (N + 1)-th ion. Therefore, we must "integrate out" these atomic degrees of freedom, which can be accomplished by projecting the density matrix of the full system intoPρ(t) = ρ ss N+1 ⊗ ρ N (t), where ρ ss N+1 fulfills D 0 (ρ ss N+1 ) = 0, and ρ N (t) is the reduced density matrix for the N spins and the chain vibrational modes. For this we use the expression where the "hats" refer to the interaction picture with respect to H 0 = H τ +H τ + H p . After some algebra, one gets the effective damping of the longitudinal modeṡ where we have introduced the bosonic dissipator together with the effective rates for laser cooling and heating κ ∓ n = D n + S n (±ω nz ).
Here, the so-called diffusion coefficient accounting for the recoil heating is D n =Γ nz τ + N+1τ − N+1 ss , and the spectral functions of the sideband terms are which can be obtained by means of the quantum regression theorem. The formalism and cooling rates coincide with those of a single trapped ion [38], with the difference that in the diffusion coefficient D n and forces F n one must consider the normal vibrational modes at the position of the cooled ion.
The idea is to work in the resolved sideband regimeΓ ω z , and set the laser-cooling parameters to optimize the cooling of the longitudinal center-of-mass mode κ − 1 κ + 1 . In this regime, the cooling dominates in (57), and we are left with the desired damping of the bosonic mode Let us finally comment on the required cooling rates. Ideally, we should achieve the regime κ ≈ g ≈ 0.1-0.5 kHz. Since such cooling is not particularly fast, it is possible to find the right detuning∆L z < 0 such that κ − 1 ≈ 0.1-0.5 kHz. Besides, since we are in the resolved sideband limit, the heating can be made much smaller, κ + 1 κ − 1 . The stationary state of a harmonic oscillator under the action of D ti b in Eq. (60) is a thermal state with mean number of quantan 1z = κ + 1 /(κ − 1 − κ + 1 ). In the limit κ + 1 /κ − 1 = ζ 1, such that the cooling is the dominant effect and the heating only presents a small correction, the equilibrium state of the center-of-mass mode is close to the vacuum state.

E. Numerical analysis of the trapped-ion dissipative protocol
In this last section, we will explore numerically how the dissipative protocol described in the part II of this manuscript can be implemented with the trapped-ion DSCB model in Eqs. (48), (51) and (60), namely Therefore, the results discussed below correspond to the actual spin-spin couplings J α i j in an ion trap, which are not nearest-neighbors couplings, but display a dipolar decay with the cube of the distance. Furthermore, we have taken into account that in harmonic ion traps, the interparticle distance is not constant over the chain. More importantly, we have considered the effects of the always-present heating term in the trapped-ion effective damping (60).
In Fig. 9(a), we show the optimal fidelity for the dissipative generation of the ideal ground state at different fillings (35), by integrating numerically the dipolar and inhomogeneous trapped-ion DSBC model (61) as a function of the ratio between heating and cooling rate ζ . We study a chain of N + 2 = 6 ions, in which N = 4 have a hyperfine structure and play the role of the spins in the XY chain, while the two peripheral ions are assumed to be used for sympathetic cooling (this choice of auxiliary cooling ions at the ends of the chain makes the interparticle spacing slightly more homogeneous for the "system" ions). As can be seen in the figure, when ζ → 0, the fidelities with the target ground state approach 100%, which allows us to conclude that the differences due to the inhomogeneity of the chain and the dipolar range of interactions with respect to the ideal homogeneous and nearest-neighbor DSBC (1) do not compromise the fidelities severely.
The situation is different as finite heating rates ζ > 0 are considered. In Fig. 9(a), we observe that the fidelity with the target ground state decays essentially linear with ζ . In the worst case considered, where ζ = 0.2 leads to a mean phonon number ofn 1z 0.25, the fidelities obtained are reduced to roughly F t f ≈ 0.8. We note that resolved sideband cooling of single vibrational modes has been experimentally achieved reachingn 1z below 0.1 for single trapped ions [36]. Sympathetic cooling of crystals up to N = 4 ions has also been achieved [26], reaching values as low asn 1z ≈ 0.06 with pulsed techniques. We note that our scheme does not demand such an accurate ground-state cooling (see Fig. 9(a)), and we expect that the required mean phonon numbersn 1z ≈ 0.1 can also be achieved in a continuous cooling scheme.
Another source of experimental imperfections is given by slow drifts of the trap frequencies which will modify the detuning δ 1z in different experimental runs, such that the conditions (19) will not be perfectly fulfilled. Nonetheless, in Fig. 4 we have shown that the dissipative character of the protocol endows it with a natural robustness with respect to nonoptimal choices of the parameters. Therefore, one can expect large fidelities even in the presence of slow drifts of the trap frequency within the kHz-range. In any case, to alleviate imperfections caused by both magnetic-field dephasing and trap frequency drifts, we have considered a faster protocol where the required time ist f ∼ 10 2 /J ≈1-10 ms. In Fig. 9(b), we show that the achieved fidelities in this case are still above 0.8 (for heating/cooling ratios below 0.2). Moreover, in this figure we have also considered substituting the standing-wave force in Eq. (51) by a traveling-wave configuration, which is experimentally less demanding.
The need of the τ z state-dependent force (49) forbids the use of magnetic-field insensitive states (i.e. clock states). Therefore, the hyperfine spins will be subject to fluctuations induced by non-shielded external magnetic fields, typically leading to dephasing in a timescale of T 2 ≈1-10 ms. These timescales are comparable to the protocol times t f ∼ 10 2 /J ≈1-10 ms, where we have taken J/2π ∼1-10 kHz. Nevertheless, these magnetic-field fluctuations act globally for standard radio-frequency traps H n = 1 2 ∆ω 0 (t) ∑ i τ z i , where ∆ω 0 (t) is a fluctuation of the resonance frequency due to the Zeeman shift (note that this might not be the case for microfabricated surface traps, where fluctuating magnetic-field gradients can also arise). Since the dynamics of the DSBC (61) conserves the number of spin excitations, this magnetic-field noise only introduces a global fluctuating phase, and thus does not decohere the state of the system. Therefore, as far as the induced spin dynamics occurs within any subspace with a conserved number of excitations, the dissipative protocol is robust to global magnetic-field noise.
It is, however, important to stabilize the intensity of the state-dependent forces to the sweet spot where J x i j = J y i j , and make sure that |F α i j | δ nα to neglect any term that does not conserve the number of excitations. The effect of off-resonant terms in the scheme leading to the XY interactions, that could give rise to non-conservation of the number of spin excitations, should be much smaller than in stroboscopic proposals [19] since our protocol does not require fast gates [45]. In Fig. 10, we show the achieved fidelities for anisotropies in J x i j /J y i j = 1. The results have been obtained numerically for the same regime as in Fig. 9(b), and indicate that if J x i j , J y i j differ by 1 part in 10 3 , the fidelities are not severely reduced. The results in Fig. 10 also show a marked difference in the robustness depending on the number of spin excitations in the initial state. This behaviour, at first surprising, can be understood by looking at the small undesired Hamiltonian terms responsible for the anisotropy in an interaction picture with respect to the ideal isotropic Hamiltonian. The anisotropy term creates or destroys fermionic quasiparticles in pairs, and rotates with a frequency that is the sum of the energies of the two quasiparticles. For the case with only one spin up, it is possible to create from the ideal target state pairs of fermions with energy adding up to zero, so that these error terms do not rotate. In the case with two spin excitations, on the contrary, the lower half of the fermionic spectrum is filled, so that any Hamiltonian term creating or destroying a pair of fermions from the ideal target state oscillates in time, and therefore its effect is strongly suppressed.

IV. CONCLUSIONS AND OUTLOOK
We have presented a method to dissipatively generate multipartite-entangled states corresponding to the ground states of small spin chains with isotropic XY Hamiltonian in a transverse field. The protocol for dissipative state preparation is interesting from a fundamental point of view, as it illustrates how local and even purely Markovian noise can assist entanglement production. Indeed, the jump operators required are sums of terms acting on individual sites, in contrast with the few-body (quasi-local) nature of the operators in [15,17]. When the noise deviates from exact Markovianity, the required time for achieving the steady state is reduced, which illustrates the value of non-Markovian effects for practical purposes. We would like to emphasize that the general idea underlying the procedure presented in this work is not specific to the model considered, and might be generalized to other spin Hamiltonians. As an example, we introduced a variation that can be used to generate states locally equivalent to W-states. It would also be interesting to modify our dissipative protocol to prepare the ground state of a gapped spin model. For longer ion chains, the combination of more state-dependent forces (4) with different wavevectors could improve the scalability of the protocol.
We have also explained in detail how to implement this method in small chains of trapped ions. In our proposal, two internal levels of the ions embody the spin system and a collective motional mode represents the damped oscillator. The preparation procedure requires the implementation of spin-spin interactions using the so-called Mølmer-Sørensen scheme, the action of a state-dependent force, and sympathetic cooling, all ingredients within the capabilities of present iontrap technology. Numerical simulations including different sources of errors indicate that the protocol can produce the target states with fidelities comparable to the more standard coherent protocols [27]. The method presented is indeed robust against a number of experimental imperfections, and the implementation is simplified since only global addressing of the ions is required.
Our results, however, may find a realization in different experimental setups. For instance, the application of our ideas to arrays of superconducting qubits in stripline resonators seems feasible. Morever, since this model can also be understood as a quadratic fermonic model with an additional chemical potential fixing the number of particles, the results may also be interesting for fermionic atoms in optical lattices.
The realization of this kind of dissipative system paves the way for the implementation of more complex scenarios to study the interplay of coherent and incoherent dynamics giving rise to noise-induced criticality [46]. Indeed, recent work in the area of non-equilibrium quantum phase transitions [47] has shown a number of fascinating results that could be demonstrated in systems of trapped ions using tools similar to the ones in our protocol.