Scalable Dissipative Preparation of Many-Body Entanglement

We present a technique for the dissipative preparation of highly entangled multiparticle states of atoms coupled to common oscillator modes. By combining local spontaneous emission with coherent couplings we engineer many-body dissipation that drives the system from an arbitrary initial state into a Greenberger-Horne-Zeilinger state. We demonstrate that using our technique, highly entangled steady states can be prepared efficiently in a time that scales polynomially with the system size. Our protocol assumes generic couplings and will thus enable the dissipative production of multiparticle entanglement in a wide range of physical systems. As an example, we demonstrate the feasibility of our scheme in state-of-the-art trapped-ion systems.

In this work, we demonstrate the range of the dissipative approach by developing general techniques for the dissipative preparation of stable and high-fidelity many-body entanglement. We show that by adding tunable sources of dissipation to a generic system of atoms coupled to a harmonic oscillator, it is possible to engineer complex manybody dissipation which can be used to generate multi-particle entangled states. We exemplify our technique considering the preparation of the most paradigmatic classes of quantum states exhibiting multi-particle entanglement, the so-called Greenberger-Horne-Zeilinger (GHZ) and W states [48], Our scheme for state preparation is deterministic and performed by continuous optical driving from an arbitrary initial state towards a desired steady state and is achieved using weak classical fields. The preparation times of our protocols are found to exhibit a favorable polynomial scaling with the number of qubits. In addition to our generic system-independent protocol we describe an implementation in a system of trapped ions. Our findings are relevant for quantum information tasks performed today in a range of physical systems which require many-body entangled states. The presented techniques are general and thus applicable to other relevant classes of quantum states and comprise a stepping stone for further advanced protocols within the framework of dissipative quantum computing and dissipative state engineering.

DISSIPATIVE MANY-BODY OPERATIONS
The challenge for any dissipative protocol for many-body state preparation is to pump an exponential number of states efficiently, i.e. in polynomial time. In the following, we demonstrate that this is possible using only a small number of elementary operations, describing the preparation of GHZ and W states: Preparation of a GHZ state (cf. Eq. (1)) starting from any initial state can, regardless of the number of qubits N , be divided into the two operations shown in Fig. 1 a): (i) Pumping all states with neither none nor all qubits in state |1 to the state |0 ⊗N , and (ii) removal of the GHZ state with the "wrong" phase, |GHZ − = (|0 ⊗N − |1 ⊗N )/ √ 2. Preparation of the W state (cf. Eq. (2)) can in principle be performed by (iii) a single collective jump from |0 ⊗N to |W (see Fig. 1 b). In order to guarantee |W to be the steady state for any initial state it is, however, necessary to also (iv) empty all states with more atoms in |1 and (v) to remove all superpositions that do not exhibit the same permutation symmetry as the W state. The mechanisms (i)-(v) are nontrivial, as they require N -body dissipative interactions. Below we present techniques to realize such operations dissipatively and thereby accomplish (i)-(v) in a realistic physical system. (a) Preparation of a GHZ state can be realized by two operations: (i) Pumping of states with 0 < n1 < N atoms in state |1 into |0 ⊗N , which is a superposition of |GHZ and |GHZ− , and (ii) removal of |GHZ− . (b) W state preparation can be accomplished by pumping all states with n1 > 1 atoms in |1 into those with n1 ≤ 1 and pumping the antisymmetric states with n1 = 1 into |0 ⊗N . Collective decay from |0 ⊗N is used to prepare |W from |0 ⊗N .

GENERIC PHYSICAL SYSTEM
For the realization of our protocols we assume a general system of N particles ("atoms"), shown in Fig. 2 (a)-(e) (see Supplementary Section I). Each atom supports two stable ground states |0 and |1 and two excited states, |e and |f . The atoms collectively couple to two harmonic oscillator modes, b and c, such as two resonator modes in cavity [14] or circuit QED [12,18,19], or two phononic modes in an ion trap setup [10,[15][16][17]. We consider two dissipative processes: Spontaneous emission from the excited states to the ground states acts incoherently on each atom. In addition, we assume loss of excitations from the oscillator modes, such as photon loss in cavity QED or phonon loss due to sympathetic cooling in a chain of ions. For our schemes we use classical fields which ideally couple identically to all atoms. In this way, we create many-body dissipation. While the setup with two excited states and two oscillator modes is completely steady, it is also possible to implement our schemes in a stroboscopic manner. Here, a single excited level and a single oscillator mode are used in two interchanging coupling situations, resulting in a quasi-steady state.

PREPARATION OF GHZ STATES
We start by describing the first dissipative mechanism (i) for the preparation of a GHZ state of N atoms. In the following, we refer to it as "Z pumping", since it is best understood in the eigenbasis of σ z = |0 0| − |1 1|, given by |0 and |1 . The process transfers all states with 0 < n 1 < N atoms in |1 towards |0 ⊗N = (|GHZ + |GHZ − )/ √ 2. We engineer the Z pumping using the coupling configuration in Fig.  2 b). Here, a coupling g of the oscillator b to the transition |e ↔ |1 and a weak drive on the same transition are used to effectively "count" the number of atoms in state |1 : As illustrated in Fig. 3 a) for the example of N = 3 qubits, the weak driving tones Ω (F ) Z couple any ground state with n 1 atoms in |1 to an atom-excited state. This state is in turn coupled to an oscillator-excited state with n 1 atoms in |1 by a coupling strength of √ n 1 g. Hence, the atom-and the oscillator-excited state form dressed states at energies ± √ n 1 g. We pump all ground states with 1 ≤ n 1 ≤ N − 1 towards those with n 1 − 1 by exciting them by weak driving tones Ω The excited states with n 1 then decay towards ground states with n 1 − 1 by spontaneous emission from |e to |0 . This operation continuously pumps population from all states with 1 ≤ n 1 ≤ N − 1 into |0 ⊗N , resulting in the production of |GHZ and |GHZ − .
To achieve the second process (ii) for the dissipative preparation of a GHZ state, we engineer the depumping of the undesired state |GHZ − by the so-called "X pumping", shown in Fig. 2 c) and Fig. 3 b). Here it is convenient to express the system state in terms of the eigenstates of σ x = |1 0|+|0 1|, |± = (|0 ± |1 )/ √ 2. Expressed in this basis, |GHZ − is a superposition of states with odd numbers of atoms in |− , n − , as written out in Fig. 3 b), while |GHZ contains only states with even n − . To depump |GHZ − we therefore need to distinguish whether n − is even or odd. We achieve this in a way similar to the Z pumping above, using the "X" coupling configuration in Fig. 2 c). Here, the transitions between the excited level |f and the ground levels |0 and |1 are coherently coupled to the harmonic oscillator c and excited by a multitone drive Ω (F ) X . Opposite phases on both transitions result in a coupling of the transition |f to |− , similar to the coupling of |e to |1 in the Z configuration. We can therefore use the X configuration to count atoms in |− in a similar manner as for |1 above. Applying fields with detunings ∆ (F ) F g for F = 1, 3, 5, . . . (F ≤ N ) thus excites the state |GHZ − and makes it decay to random states by spontaneous emission, as shown in Fig. 3 b). On the other hand, |GHZ remains unaffected. After decay out of |GHZ − , the resulting states are again subjected to the Z pumping (i) and thus transferred back to |0 ⊗N . The combination of both processes thus dissipatively prepares |GHZ over time and maintains it as the unique steady state of the dissipative dynamics (see also Supplementary Section IV). However, since the X pumping disturbs the Z pumping, it has to be sufficiently weak so that the Z pumping has a sizeable probability of reaching the final state |0 ⊗N before being subject to X pumping; this requirement does not, however, slow down the preparation process significantly (see Supplementary Section VII). a) FIG. 2. System. (a) We consider a chain of N sub-systems ("atoms") with four levels, coupled to two common harmonic oscillators. As dissipative processes we assume spontaneous emission from the excited levels γ e/f and oscillator decay κ b/c . (b)-(e) Depending on the scheme, we apply coupling configurations ('Z', 'X', 'A', 'W') consisting of classical drives Ω with detunings ∆ and couplings of the atoms to the oscillator modes g (see also Supplementary Section II).

Scaling of the protocol
The resulting time-evolution towards a steady GHZ state for qubit numbers N = 2, . . . , 8 is shown in Fig. 4 a). The curves are obtained by numerically simulating an effective master equation (see Supplementary Section III), after optimization of the available system parameters to achieve a desired fidelity of F GHZ = 0.9 in minimal time. The inset shows the time required to reach this fidelity when starting from an initially fully mixed state, as a function of the system size. It can be seen that the preparation time behaves better than a polynomial bound which we derive from a simplified analytical model (see Supplementary Section VII): In short, we find that the steady state fidelity F is determined by the ratio of the decay out of the GHZ state at a rate Γ − ∼ γN Ω 2 /g 2 , and the effective preparation rate Γ + ∼ Ω 2 /γ, giving a steady state error 1 − F = Γ − /Γ + ∼ N γ 2 /g 2 . With a suitably chosen low decay rate γ ∼ g √ E/ √ N we can prepare a GHZ state with any desired error 1 − F = E. Such an adjustable decay rate can, e.g., be achieved by taking |e and |f to be metastable states in an atom and driving them to an excited state. By adjusting the driving we can thus vary the decay rate [45,46]. The preparation rate Γ + increases with the strength of the driving Ω, but for too strong driving the fidelity becomes affected by power broadening, limiting the preparation rate (see Supplementary Section VI). As shown in Supplementary Section VII, these detrimental effects can be suppressed by a suitably low driving strength Ω ∼ γ/ √ N . In this way, we derive a rigorous upper bound on the preparation time whose scaling behaviour is given by Our bound is shown in the inset of Fig. 4 a), along with the numerical results.

PREPARATION OF W STATES
To realize the protocol in Fig. 1 b) for the preparation of W states we again use the Z pumping in Fig. 2 b) and Fig. 3 a), together with two more coupling configurations, "A" and "W", shown in Fig. 2 d) and e). For the Z pumping we now apply fields with detunings ∆ are pumped to |0 ⊗N by the "A pumping": As shown in Fig. 2 d), we apply a resonant drive (∆ = 0) to the transition from |1 to |f and an oscillator coupling on the transition from |f to |0 . Ground states of the n 1 = 1 manifold are then driven to atom-excited states of the same symmetry. The fully symmetric atom-excited state to which |W is excited by the drive, couples to the state |0 ⊗N |1 c which has an excitation of oscillator c. Due to their coupling of √ N g these states form dressed states at energies ± √ N g, as illustrated in 3 c). Excited states of other symmetries, to which the antisymmetric states are coupled, are dark states of the oscillator coupling and remain in resonance with the drive. The antisymmetric ground states are thus excited rapidly to the corresponding atom-excited states from where they decay to |0 ⊗N by spontaneous emission. |W , in turn, remains unaffected due to the detuning of the corresponding dressed states.
The above mechanisms pump all states except for |W to |0 ⊗N . The actual preparation of |W from |0 ⊗N is then realized by the "W pumping": This mechanism is achieved by an effective decay (or "cooling") of the harmonic oscillator, realized by the coupling configuration in Fig. 2 e). Here we apply a two-tone drive on the transition from |0 to |e and use an oscillator coupling on the transition from |e to |1 . As shown in Fig. 3 d), |0 ⊗N is then excited by the drive to a symmetric atom-excited superposition state, from where the atomic excitation is transferred to the oscillator, yielding the state |W |1 . As the latter two states form dressed states at energies of ±g, choosing detunings of ∆ W = ±g ensures that |0 ⊗N is excited resonantly. Oscillator decay then transfers the system to the |W state. Due to the specific detuning, leakage of population from |W to higher Dicke states does not occur, as this would require detunings ∆ = ± √ n 1 g with n 1 > 1. The simultaneous action of the three mechanisms Z, A, and W yields the desired W state as the unique steady state (Supplementary Section V). Alternatively, the W mechanism can also be realized by spontaneous emission at the cost of a factor ∼ √ N in the preparation time.   Fig. 2 (b). Ground states are coupled to atom-excited states by weak driving. Depending on the number n1 of atoms in |1 , the excited states form dressed states with oscillator-excited states at energies ± √ n1g. By applying fields with detunings ∆ Depumping of |GHZ− ("X pumping") uses the couplings in Fig. 2 (c). |GHZ− is a superposition of states with only odd numbers n− of atoms in |− , whereas n− is even for |GHZ . The atom-excited states to which |GHZ− is coupled form dressed states at ± √ n−g with odd n−. |GHZ− is then depumped by applying fields with detunings ∆ (F ) To pump all states with n1 ≥ 2 to n1 ≤ 1 we apply the mechanism in (a) with detunings ∆ The antisymmetric states |1 k (1 ≤ k ≤ N − 1) with n1 = 1 are emptied by the "A pumping", using the coupling situation shown in Fig. 2 (d): While the atom-excited state to which |W is excited forms dressed states with |0 ⊗N |1 at ± √ N g, the atom-excited states to which the antisymmetric states with n1 = 1 are excited are dark states of the oscillator coupling. They are therefore in resonance with the drive at the bare atomic resonance and pumped to |0 ⊗N . (d) |W is prepared by "W pumping", using the coupling situation shown in Fig. 2 (e). An effective decay process consisting of weak excitation with a detuning ∆W = ±g, atom-oscillator coupling g, and oscillator decay κ leads from |0 ⊗N to |W , whereas excitation out of |W is suppressed as this would require a detuning of √ 2g.

Scaling of the protocol
In Fig. 4 b) we give the scaling of the preparation time of the W protocol found by numerically simulating the effective dynamics of the scheme (see Supplementary Section III) using the collective W mechanism and optimizing the available parameters. Using similar arguments as above (see Supplementary Section VIII), we find that for a desired error E the preparation time scales as The scheme thus yields a low-order polynomial scaling of the preparation time with the number of qubits. Note that, while here we have assumed that we initially start from a fully mixed state, starting from the state |0 ⊗N would improve the scaling of the time for W preparation significantly.

PHYSICAL IMPLEMENTATION
All the ingredients necessary for our scheme are available in trapped ion experiments: A suitable setup can consist of a chain of N trapped ions, with two (meta-)stable ground levels |0 and |1 and two auxiliary levels, |e and |f . Two phononic modes, cooled to the ground state, e.g. by sympathetic cooling of auxiliary ions, and coupled to the sidebands of the ions' transitions, can be used as harmonic oscillators. Tunable decay of the auxiliary levels by spontaneous emission can be realized by a repumper to a rapidly decaying state [45,46]. An alternative stroboscopic implementation requires only a single-auxiliary level with a repumper, interchanging between the role of |e and |f in the assumed coupling configurations, and a single phononic mode, interchanging between b and c. DISCUSSION We have shown that dissipative state preparation can be extended to the generation of many-body entanglement. We  FIG. 4. Evolution towards steady many-body GHZ and W states. Starting from an initially fully mixed state, we numerically solve an effective master equation of the protocol for dissipative preparation of (a) GHZ and (b) W states. The curves show the evolution for two to eight qubits (different colors) and are obtained by numerically optimizing all available parameters to reach a fidelity of F = 0.9 of the desired state within as short as possible a preparation time. The insets show the scaling of the preparation time with the number of qubits. In the curves and in the insets, we show different degrees of truncation of the Hilbert space (dash-dots/blue squares -effective dynamics after adiabatic elimination, solid lines/green circles -one excitation, dashes/red triangles -two excitations). In the insets, small symbols stand for analytically optimized parameters and large symbols for numerically optimized parameters. We find favorable polynomial scalings both for the GHZ and W preparation times which are within our analytical bounds (black dashed lines).
achieve this by engineering multi-particle dissipation which deterministically drives the system into a desired steady state within a time scaling only polynomially with the size of the system. Our approach may be implemented in physical systems such as trapped ions where the basic ingredients of the scheme have already been demonstrated [46].
We have chosen to exemplify the procedures by preparing GHZ and W states, which are the paradigmatic examples of multipartite entangled states, but the developed techniques are applicable to a range of other tasks. Particularly interesting possibilities are the construction of quantum error correcting codes [2][3][4]50] or the observation of exotic phase transitions [24] induced by multi-particle dissipation.
The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007(FP/ -2013

SUPPLEMENTARY INFORMATION
In this Supplementary Information to the article "Scalable dissipative preparation of many-body entanglement" we present the dynamical model of the generic system under consideration (Section I), the coupling configurations for the described protocols (Section II) and the resulting effective dynamics of the system (Section III). We discuss the engineering of the dissipative manybody interactions for GHZ and W state preparation (Section IV and Section V) and strong driving effects (Section VI). In Section VII and VIII we analyze the scaling of the protocols both for weak and for strong driving.
Here, ρ is the density matrix of the system and L denotes a Liouvillian of Lindblad form. The Hamiltonian H contains the coherent interactions, while sources of dissipation present in the system are represented as Lindblad ('jump') operators L k . The Hamiltonian for the coupling configurations detailed in Fig. 2 b)-e) is, in its most general form, given by We consider a free Hamiltonian H free which contains the energies of the levels of the N atoms and the harmonic oscillator modes, Here we have introduced J ij = N a=1 σ ij a = N a=1 |i a j| and made the simplifications ω 0 = ω 1 = 0 and = 1. An interaction Hamiltonian describes the atom-oscillator coupling, and a drive Hamiltonian H drive contains the fields used to perform coherent excitations of the system. The interaction terms and drives required for GHZ and W preparation are detailed in Section II.
The excited degrees of freedom in the system are generally subject to dissipation. Here the excited states |e and |f of each atom undergo spontaneous emission to each of the ground states |0 and |1 , described by the jump operators where the subscript a denotes the atom number. The total decay rates of the excited levels are given by γ e = γ 0e + γ 1e and γ f = γ 0f + γ 1f . For simplicity we will assume equal decay rates to both ground states, which is, however, not crucial for the protocol. The decay of excitations of the two oscillator modes, b and c, is represented by The GHZ scheme does not require oscillator decay at all (κ b = κ c = 0). It may, however, still be useful to avoid heating. For W preparation, on the other hand, we assume dissipation of one of the modes (κ b > 0).

II. COUPLING CONFIGURATIONS
The atom-oscillator couplings of the four coupling configurations in Fig. 2 b)-e) of the main text are described by the interaction Hamiltonians Here, 'Z', 'X', 'W' and 'A' denote the coupling configurations (recall that |− = 1 √ 2 (|0 − |1 )). For GHZ preparation we use the Z and X configuration, i.e. H int = H int,Z/W + H int,X , and for W preparation the Z, W and A configurations, i.e. H int = H int,Z/W + H int,A . These terms describe that an atomic excitation can be exchanged coherently with the respective harmonic oscillator with a coupling constant g. By H int,Z/W an atomic excitation in |e is exchanged with the oscillator b, leaving the atom in |1 . H int,X couples the excited level |f to |− = 1 √ 2 (|0 + |1 ) while exchanging the atomic excitation with the oscillator c. The coherent excitation of the atoms by classical driving fields is modeled by the drive Hamiltonians Here, we generally allow for several field tones with Rabi frequencies Ω (F ) k and frequencies ω (F ) k , from which we will deduce detunings ∆ (F ) k of the respective fields. For GHZ preparation we use the Z and X configuration, i.e. H int = H int,Z/W + H int,X , and for W preparation the Z, W and A configurations, i.e. H int = H int,Z/W + H int,A .

III. EFFECTIVE DYNAMICS OF THE OPEN SYSTEM
In the following, we provide a detailed analysis of the effective dynamics of the system described in Section I. To this end, we briefly introduce an effective formalism presented in Ref.
[49] and use it to derive effective operators for the coupling configurations presented in Section II.

A. Effective operator formalism
As can be seen from Section I, the dissipation affects the excited levels |e and |f and the oscillator modes a and b. For weak driving the decaying degrees of freedom can be adiabatically eliminated from the master equation. This is done using the effective operator formalism presented in Ref. [49]. In this way, the dynamics of the master equation are reduced to effective couplings between the ground states of the system, described by an effective master equatioṅ Since we are dealing with multiple field tones F that give rise to the effective couplings, we use the extended formalism for many fields (cf. [49]) with The non-Hermitian Hamiltonian H (F ) NH , which contains the frequency ω (F ) , describes the time evolution of the excited states. The drives are regarded as perturbations and thus denoted by "V "; they are defined by the drive Hamiltonians H drive in Section II so that we use V = H drive when we derive the effective operators for each of the coupling configurations below. Based on the assumption of weak excitation we will restrict our discussion to the states which have at most one atomic or oscillator excitation. H NH then contains the energies and couplings of the excited states. Since H NH needs to be inverted to compute the effective operators, we will, for each coupling situation, start out by discussing this entity.
B. Derivation of the effective operators for the coupling configurations

Z configuration
We begin with the Z coupling configuration, which is similar to the X configuration. Both are required for the generation of GHZ states. We use H int,Z , H drive,Z from Eqs. (S14)-(S15) and the Lindblad operators in Eqs. (S8)-(S13) to set up the non-Hermitian Hamiltonian Here we have introduced complex detunings∆ where F denotes the particular tone of the driving field. Furthermore we have changed into a frame rotating with the frequency of the drive ω (F ) Z . For the derivation of the effective operators we will for simplicity drop the sub-and superscripts denoting the coupling configurations and field tones. To invert the non-Hermitian Hamiltonian we divide it into four blocks After this separation we can formally invert the Hamiltonian H NH , using Banachiewicz' theorem [51] for the blockwise inversion of a square matrix, We now need to compute d to obtain any of the above elements. As we shall see, it is possible to simplify the calculation and to obtain closed expressions for the decay rates if we separate the involved operators by the number of atoms in state |1 . This is done by introducing projection operators P n1 which project on the states with the same number n 1 (from now on, n) of atoms in |1 . For the Z and X configurations discussed here, n 1 (or n − for X), is conserved under the couplings by the coherent interactions H int and V , but can be changed by the dissipative jump processes L k , e.g. from n 1 to n 1 − 1 in the case of Z. Using these projectors we can split the non-Hermitian Hamiltonian of the excited states and its four blocks by n, The inverse of these non-Hermitian Hamiltonian for each n is then found to be The effective operators of Eqs. (S22)-(S23) are formally given by To obtain the effective Lindblad operators and the effective Hamiltonian it is thus sufficient to compute the blocks d n and b n of the inverse non-Hermitian Hamiltonian. Using the identities a(a † a) −1 a † = 1 and P g J 1e J ee = P g J 1e (where P g is the projector onto the ground states) we obtain With this, and readopting the sub-and superscripts for the configuration and the field, and changing to a more detailed notation for the effective Lindblad operators we find To obtain this we have used the identities P g J 1e J e1 P g = P g J 11 P g , J 11 P n1 = n 1 P n1 , and J ee J e1 P g = J e1 P g . In the last expression we have also introduced the effective detunings∆ Z,n1 and the effective couplingsg Z,n1 with As can be seen from Eqs. (S38)-(S40), dealing with multiple frequencies in the drive leads a priori to a sum over terms for all fields in the effective Lindblad operators. However, as the frequencies of these fields are well-distinguishable, we separate the Lindblad operators by their driving field F . Given the quadratic appearance of the Lindblad operators in the master equation, we can also drop the exponential phase factors. For the effective Lindblad operators for the fields F we then obtain We also define the corresponding effective decay rates The operators in Eqs. (S43)-(S45) are then the effective Lindblad operators for the Z configuration. As can be seen from the expressions in Eqs. (S41)-(S42), the effective detunings∆ (F ) Z,n1 can be engineered to be very small by a suitable choice of the frequencies ω (F ) Z of the fields F which can be used to engineer the rates γ (F ) 0,Z,n1 and γ (F ) 1,Z,n1 of the effective decay processes. The engineering of the effective decay process to prepare GHZ and W states will be subject to Section IV. We now turn to the effective Hamiltonian. The effective Hamiltonian is computed from Eq. (S22). Other than for the effective Lindblad operators, introducing a multi-tone driving field results in cross terms between different fields, here denoted by F and G, Here, all terms F = G have fast rotating exponential phase factors. Restricting the treatment to F = G where these terms cancel, we obtain the main contribution We thus find that the main effective Hamiltonian processes are AC Stark shifts with a magnitude As we will see further below, our choice of the field tones will make these Hamiltonian terms compensate each other.

X configuration
We perform an analogous treatment for the X pumping meditated by the excited level |f and the oscillator mode c. The effective operators are obtained in the same way as for the Z configuration, using the non-Hermitian Hamiltonian Carrying out the derivation in the same manner as above for Z, we obtain for the effective Lindblad operators with the effective detunings∆ The effective decay rates can be written as These operators resemble the ones for Z pumping in Eqs. (S43)-(S45) if |1 is replaced by |− . The only difference is that spontaneous emission is still assumed to lead to the final states |0 and |1 .

W configuration
We now turn to the other two coupling situations, W and A, which are required for preparing the |W state. The situation for W and A is different from the previous cases of Z and X in the sense that H int and H drive act on different transitions. Here, we first consider W and, in the next section, A. To find d (F ) n we again consider H NH in blocks, denoted by the number of atoms in state |0 , n 0 , that can be excited by the drive. Starting from and using∆ Here, we have defined the effective couplings With the identity P g J 1e J e0 P g = P g J 10 P g this yields the effective operators for W, It can be seen that for the effective operators for the W configuration -and, as we will see below, also for A -we typically do not reach expressions which are as compact as for the Z and X configuration. However, as we will see below, these expressions simplify considerably for most of the states which are relevant to the scheme, in particular for the dark states of J 10 .

A configuration
The effective operators for A are found by making the replacements 0 ↔ 1 and e ↔ f and n 0 ↔ n 1 in the terms for W above.
We then obtain In Sections IV and V we will engineer the obtained effective operators by a suitable choice of the parameters.

C. Reduction of the dynamics to rate equations
Using the effective operator concept, we have so far reduced the dynamics of the open system to an effective master equation of the ground states. This description is exact to second order perturbation theory. In addition to that, we will later achieve a compensation of the effective Hamiltonian, H eff = 0, so that the remaining dynamics are purely dissipative. We can then, in another step, reduce the complexity of the dynamics to rate equations of the populations. This is achieved by choosing subspaces of interest between which the interactions present in the system do not build up coherences. These subspaces are defined by projection operators, e.g. P A and P B . For negligible coherences between the subspaces, we can then trace over the Liouvillian evolution to obtain decay rates from subspace A to subspace B For the subspaces we will consider, the decay rate is the same for all states in the subspace. We can then calculate the decay rates using a single state |ψ i , We will use this rate equation approach to analyze the scaling of the preparation time of the protocols in Section VII.

D. Numerical simulations
The evolution of the system is simulated by numerically solving the effective master equation (S21), including the effective operators derived above. For the evolution due to the effective master equation of the GHZ protocol we use a Trotter-like ansatz, simulating the evolution under the Z and X coupling in an interchanging manner, performing base transformations between the eigenbases of σ z and σ x in-between. The simulations for the W scheme are carried out entirely in the σ z basis. To verify our findings, we compare our result to the solution of the master equation (S5) after truncation to one or two excitations.

IV. ENGINEERING DISSIPATIVE MECHANISMS FOR GHZ STATE PREPARATION
In this section we show how the effective operators derived in the previous section are engineered to prepare GHZ states efficiently. As is discussed in the main manuscript, the engineered operators are used to empty certain states and transfer the population to others in such a way that in the end only the target state remains. This is achieved using the driving fields in the various coupling configurations Z and X (and Z, W and A in case of W state preparation) which activate the effective decay processes derived in Section III. By choosing suitable detunings for these driving fields we can then engineer the rates of the effective decay processes. The choice of the Rabi frequencies of these fields is subject to Section VII where we optimize the preparation time of the protocols.
A. Preparing |GHZ : the Z pumping The first mechanism to prepare a |GHZ state described in the main manuscript is the Z pumping: This process transfers the population of all states other than |0 ⊗N and |1 ⊗N to |0 ⊗N = 1 √ 2 (|GHZ +|GHZ − ), and thereby, to |GHZ , without affecting the GHZ state. We engineer this process by applying drives Ω Here, the subscript Z+ denotes the fields with positive detunings and Z− those with negative detunings. If the field index F coincides with the number of atoms in |1 , n 1 , for a certain initial state and driving field, the field is resonant and the effective detunings in Eqs.
Since we generally work in the strong coupling limit γ, κ g, the above effective detunings are small compared to those for off-resonant driving fields with n 1 = F ,∆ With the effective detunings we can then compute the effective decay rates in Eq. (S46) -(S48) for the processes with F = n 1 and for the off-resonant ones with F = n 1 , From these expressions we can clearly see that the first group of rates is engineered to be strong, while the second group of rates is engineered to be suppressed. Using the quantities above we obtain for the effective Lindblad operators in Eqs. (S43)-(S45).
From these expressions we see that for 1 ≤ n 1 ≤ N − 1 we will always find terms with n 1 = F to zeroth order in g −1 , which are much larger rates than the ones with n 1 = F , which are to second order in g −1 . We can therefore drop the latter for 1 ≤ n 1 ≤ N − 1. The terms with n 1 = N which in particular affect the GHZ state need, however, to be kept. Since the effective detunings are engineered to depend on n 1 , photons scattered by resonances with different n 1 can be distinguished. Formally this is justified by the exponential factors e −iωt , washing out interferences between terms (cf. Section III). We can therefore separate the terms with different n 1 into individual Lindblad operators, each acting on a set of states with n 1 atoms in |1 . The enhanced processes are then given by The weak decay processes affecting |1 ⊗N , and thus, |GHZ are found to be The result of the engineering of the effective processes for the Z configuration is thus the Z pumping by which all population from states with 1 ≤ n 1 ≤ N − 1 is transferred to |0 ⊗N , and thus, to |GHZ ± . Beside effective decay processes, we generally also obtain effective Hamiltonian terms as given by Eq. (S51), which we write as These terms are AC Stark shifts, the magnitude of which is obtained from Eq. (S52), By our combination of red-detuned drives (Ω (F ) Z− we achieve that the shifts compensate each other, Given that there are no Hamiltonian terms, we can turn to a description of the dynamics in terms of rate equations: In Section III C, we have described the possibility to reduce the effective dynamics further to rate equations. To perform such a step, we identify subspaces inside which all states have the same decay rate and in-between which no significant correlations are built up. It is therefore important that the effective Hamiltonian is zero. We define the subspaces by the projectors P GHZ = |GHZ GHZ|, P GHZ− = |GHZ − GHZ − |, and P n1 , which contain the states with a certain number of atoms in |1 . The decay rates can then be calculated using Eq. (S77). For resonant Z pumping (F = n 1 ) from a subspace with n 1 atoms in |1 to one with n 1 − 1 due to spontaneous emission to |0 we obtain the rate The subscripts in the above and the following decay rates specify the initial subspace, the final subspace, the physical process, i.e. oscillator decay (κ), spontaneous emission to |0 (γ 0 ) or |1 (γ 1 ), and the pumping process (here: Z). As opposed to the resonant rate (F = n 1 ) above, the decay rates due to off-resonant fields can be written as (S102) From these expressions follow the loss rates from |GHZ due to Z pumping: (S105) Here we have neglected the gain of population in |GHZ from |GHZ − . Furhtermore, we note that due to its scaling with N 2 , loss from |GHZ by oscillator decay should be avoided. In addition, effective oscillator decay is not useful for the Z pumping process. As we will see below, this is also the case for X pumping. Therefore, we will generally choose to work with κ b/c = 0 for the GHZ protocol (a weak cooling κ Ω may nevertheless be used to counteract heating). The overall loss rate from |GHZ due to off-resonant Z pumping is then given by Here, the question mark stands for any potential final state; in the scaling analysis we will typically consider the worst state possible. For the reasonable assumption of γ 0e = γ 1e = γe 2 we obtain We conclude that for g γ, κ, the rates from Z pumping, ultimately leading to |0 ⊗N , and thus to |GHZ , are much stronger than the loss rates from |GHZ . The derived rates will be used further to analyze the error and preparation time of the protocol in Section VII.
B. Emptying |GHZ− : the X configuration |GHZ − is emptied by X pumping, as described in the main part: Here we make use of the fact that in the eigenbasis of σ x = (|1 0| + |0 1|), consisting of |± = 1 To depump |GHZ − without affecting |GHZ we apply again red-and blue-detuned laser fields Ω For the effective decay rates we find from Eqs. (S59)-(S61), for resonant excitation (F = n − ) and for off-resonant (F = n − ) excitation. With these rates and Eqs. (S54)-(S56) we obtain the effective Lindblad operators 1,X±,n− |1 a −|P n− , (F = 1, 3, 5, ..., (F ≤ N )). (S117) Again, we separate the effective Lindblad operators by the frequencies of the resonances, this time depending on n − . We obtain similar Lindblad operators as for Z pumping, with enhanced terms to zeroth order in g, X±,n− P n− , (F = 1, 3, 5, ..., (F ≤ N )), and suppressed terms acting on the target state It can be seen that for g γ, κ these operators make |GHZ − (with only odd n − ) decay rapidly. The losses from |GHZ (with only even n − ) are caused by the off-resonant drives with F = n − . Since these terms are of second order in g −2 , they are suppressed for g γ, κ. Again, the effective Hamiltonian is compensated by the red-and blue-detuned fields, We can therefore describe the dynamics in terms of rates: With the effective operators from Eqs. (S118)-(S120) we find the decay rates from |GHZ − , (S126) For the reasonable assumption of γ 0f = γ 1f = γ f 2 the total rate is therefore given by The expression in Eq. (S127) is the total decay rate from the |GHZ − state which is approximately given by the sum of all enhanced decay rates, weighted with the number of states with the same excitation. The X pumping distributes the population of |GHZ − over all other states which we here denote by a question mark and later on replace by a worst-case assumption.
The losses from |GHZ are only caused by the off-resonant drives with F = n − . Using Eqs. (S121)-(S123) and ignoring negligible gain processes we obtain the loss rates Here, the binomial coefficients originate from the number of states with the same number of atoms in |− . For the total loss rate from |GHZ through X pumping with both red-(X+) and blue-detuned (X−) fields we approximately find From this expression we see that the loss terms from |GHZ due to X pumping are of second order in g −1 and thus suppressed for g γ, κ. Beside decay out of |GHZ − the X pumping also causes losses from states with odd n − which have an overlap with states with 1 ≤ n 1 ≤ N − 1 in the Z basis. This affects the transport from n 1 = N − 1 to n 1 = 0 by the Z pumping and thus to |GHZ by imposing a loss rate We will refer to this process as "X toss" below. As we will show below this does not, however, have a significant effect on the scaling of the preparation time and the error, if the strength of the X process is chosen properly. We conclude that the combination of Z and X pumping results in the preparation of a GHZ state from any initial state. For g γ, κ the gain rates are engineered to be strong, while the loss rates from |GHZ are suppressed. The scaling of the error and the preparation time of the protocol are analyzed in Section VII.

V. ENGINEERING DISSIPATIVE MECHANISMS FOR W STATE PREPARATION
We analyze the effective operators and effective transfer rates for the W preparation scheme. Here, we first discuss the Z mechanism used to pump states with n 1 ≥ 2 atoms in state |1 to those states with n 1 ≤ 1. Then we derive the rates for the depumping of the N − 1 antisymmetric states within the subspace with n 1 = 1. Finally, we analyze the actual preparation mechanism, the transfer from |0 ⊗N to |W . Since the effective operators for spontaneous emission in the W and A configuration Eqs. (S67)-(S69) and (S71)-(S73) do not have a simple form, we restrict our discussion of the effective operators to the manifolds with n 1 ≤ 2, i.e. states with at most two atoms in state |1 . This is justified as the W and A operations have a negligible effect on states with n ≥ 2.
A. Transferring states to n1 ≤ 1 by the Z pumping We analyze the Z pumping from states with n 1 ≥ 2 into those with n 1 ≤ 1. To do so we insert the detunings ∆ (F ) Z± = ± √ F g for 2 ≤ F ≤ N into the expressions for the effective detunings in Eqs. (S41)-(S42) and decay rates (S46)-(S48) This allows us to compute the Lindblad operators in Eqs. (S43)-(S45) for this parameter choice. For the enhanced processes (n 1 = F ) we find In addition, there are weak decay processes affecting |W described by The effective Hamiltonian is again compensated by applying drives with both blue and red detunings (cf. Eq. (S99)), so that we can turn to a description using rate equations.
We first calculate the rates for the desired mechanism consisting of several enhanced processes and then the loss rates. The decay rate leading from n 1 = 2 into n 1 = 1 constitutes an upper bound for the decay rate from n 1 = 2 to the antisymmetric subspace (from now on denoted by "as"), which contains all states |1 k =0 with n 1 = 1 except for the W state. This rate is given by (S138) Generally, the rate for Z pumping from a subspace with n 1 atoms in |1 to one with n 1 − 1 is found to be These rates will be needed to derive the time it takes to reach |W by several decay processes. We now turn to the effective decay processes from |W . Here, we find for the losses due to Z pumping As can be seen from Eq. (S139), the decay rates of the Z process depend on the number of atoms in |1 of a state, and not on its symmetry. The antisymmetric states are emptied by the A pumping discussed below.

B. Emptying antisymmetric states: the A pumping
Since the Z pumping does not distinguish between the symmetry of the states in n 1 = 1, the A pumping is used to empty all states other than |W . Here, it is helpful to recognize that, while |W is completely symmetric under permutation, the other N − 1 states with n 1 = 1 are antisymmetric under permutation. A basis for the antisymmetric states is given by where σ + a = |1 a 0|. Note that |1 0 = |1 N = |W is the W state, or first Dicke state. The second Dicke state is the symmetric superposition of all states with two atoms in |1 . The A mechanism is meant to pump the antisymmetric states to |0 ⊗N , while ideally leaving |W unaffected. To achieve this we use the A configuration with a resonant single-tone field, ∆ A = δ A = 0. The decay of the oscillator mode c is assumed to be negligible, κ c = 0. With this we obtain for the Lindblad operators Eqs. (S71)-(S73) in the manifold n 1 ≤ 1 Here we have defined a projector P n1=1,as = N −1 k=1 |1 k 1 k | onto all antisymmetric states in n 1 = 1 and a projector onto the W state, P W = |W W|, and the effective decay rates with the effective terms from Eqs. (S65)-(S66),g It can be seen that the effective spontaneous emission from antisymmetric states with n 1 = 1 is only limited by the line width γ f and thus very rapid compared to the loss terms suppressed by g −2 . The losses from |W are, in turn suppressed by a detuning growing with the qubit number and switched off completely for κ c = 0. The effective Hamiltonian of the A mechanism acting on the n 1 = 1 manifold is zero due to the absence of a detuning of the driving field, H A ≈ 0 (S151) so that we can describe the A pumping by rates. The depumping rate for antisymmetric states ("as") is given by The loss rates from |W are, on the other hand, found to be small, The A pumping is found to be limited only by the more or less freely adjustable decay rate γ f of the excited level |f . The depumping of the antisymmetric states to |0 ⊗N is thus a very efficient process. Preparation of the |W state from |0 ⊗N is performed by the W pumping discussed in the following section.
C. Preparation of |W from |0 ⊗N : the W pumping For the W mechanism we use a two-tone field with a single red and a single blue detuning, ∆ W ± = ±g = δ W± . By this condition we transfer population from |0 ⊗N (n 1 = 0) rapidly and directly to |W (n 1 = 1) by collective decay, while loss from |W is effectively suppressed by different resonance conditions. Preparation of antisymmetric states |1 k =0 is also avoided unless spontaneous emission dominates over oscillator decay, γ e κ b , in which case a modified scheme for W preparation with a different scaling is obtained. We find for the effective operators of Eqs. (S67)-(S69) in the manifold n 1 ≤ 2, Here, we have introduced the effective decay rates with effective couplings for the drive acting resonantly on |0 ⊗N (F = 1, 0), or off-resonantly on |W (F = 1, W ), Again, the effective Hamiltonian processes are compensated by the use of fields with opposite detunings, so that we can restrict the discussion to gain and loss rates. Using the engineered operators above we find for the preparation rate of |W , In addition, there are enhanced decay processes from |0 ⊗N to states with n 1 = 1 by spontaneous emission The first of these processes can be used to build an alternative W preparation scheme which works in the absence of oscillator decay, κ b = 0. Here, population pumped to "as" by the rate above is transferred back to |0 ⊗N by the A mechanism. The less directed preparation mechanism results, however, in a different scaling. For the loss rates from |W due to the W process we We find that all loss rates from the target state |W scale with 1/g 2 , while the gain terms are only limited by the decay rates γ and κ. We therefore conclude from our analytical discussion of the effective operators of the W scheme that the combination of the three presented mechanisms allows for the dissipative preparation of a W state. Since all states in the system eventually end up in the W state, |W is attained as the unique steady state of the dynamics, independent of the initial state. Since for g γ, κ the gain rates are much larger than the loss rates from |W the W state can be prepared within short time and to a high fidelity. We investigate these questions further in Section VIII, where we use the derived rates to analyze the scaling of the error and the speed of the protocol.

VI. STRONG DRIVING EFFECTS
The operator formalism used to derive effective couplings of the ground states in the system is built on the assumption of a perturbative drive which is much smaller than the couplings, decay, or detunings of the excited states, Ω γ, κ, g, δ, ∆. Effects from saturating the excited states that could possibly slow down the preparation process are thus not included in the rates derived so far. On the contrary, the dependency of all decay rates on the drives ∝ Ω 2 suggests that the optimal Rabi frequencies Ω of the drives are infinitely strong and thereby outside the perturbative regime. A proper assessment of the driving strength Ω therefore requires the inclusion of strong driving effects which begin to play a role for Ω γ, κ. In the following, we include such effects in the dynamics in an approximate way. Here, we discuss power broadening and population of the excited states and adjust the rates derived from the effective operators to account for power broadening of the excited states. They will later be used to analytically derive the optimal Rabi frequencies Ω and the scaling of the protocols for strong driving. In addition, in our numerical simulations we use effective operators where we have included the power broadening terms manually. While this treatment is not rigorous, it provides a convenient tool for rapidly determining the optimal parameters. These are used to simulate the effective dynamics of the system beyond the regime of weak driving and to match it with the simulations describing a larger Hilbert space including excitations. In addition, we take into account the effect of population of the excited states in our numerics.

A. Power broadening
We first address the effect of power broadening (or 'line' broadening), regarding a simple model situation: A ground state |0 is resonantly coupled by a field with a Rabi frequency Ω to an excited level |e . In total, |e decays at a rate γ = γ 0 + γ 1 , where γ 1 is the decay rate into |1 . We perform adiabatic elimination by setting the derivative of the density matrix to zero,ρ ≈ 0. In the weak driving regime, where the broadening of the excited level can be neglected, this yields an excited population of ρ ee = Ω 2 γ 2 ρ 00 and thus an effective decay rate from |0 into |1 of γ eff = Ω 2 2γ . The population gain of state |1 is then given bẏ The same result is obtained when using the effective operators. Performing adiabatic elimination with a stronger drive we need to take into account the population of the excited level. This yields an excited population of ρ ee = Ω 2 γ 2 +Ω 2 . Then, as the coupling of the ground state |0 to |e is increased, we need to take into account the population of the coupled subspace of |0 and |e rather than of state |0 only. This leads tȯ ρ 11 ≈ γ 1 ρ 00 (ρ 00 + ρ ee ) (ρ 00 + ρ ee ) ≈ γ 1 Ω 2 γ 2 + 2Ω 2 (ρ 00 + ρ ee ) (S172) In the weak driving limit Ω → 0, this new decay rate γ eff = γΩ 2 2(γ 2 +2Ω 2 ) matches the previous result. The decay rate for strong driving can thus obtained from the rate which was derived from the effective operators by "broadening" the natural line width of the excited state, γ 2 → γ 2 + 2Ω 2 . In the limit of strong driving Ω → ∞, the population is found in |e with a probability of 1 2 such that here γ eff approaches the constant value of γ1 2 . Therefore, the effective decay rate from |0 to |1 cannot be increased beyond half of the line width of the level mediating it. More generally, we now seek to include power broadening in the effective operators. This is done by replacing the complex effective detunings, e.g.∆ (F ) Z,n1 , here generally denoted as∆ eff andg eff , by "power broadened" ones. The replacements are made such that the rates obtained from the effective operators agree with the ones derived by adiabatic elimination. To this end, we make the replacements of |∆ eff | 2 → |∆ eff | 2 + nΩ 2 and |g eff | 2 → |g eff | 2 + nΩ 2 , where n is the number of atoms that can be excited by the drive. While our numerical simulations show that this replacement is sufficient to account for power broadening in the A operators, it turns out that the action of power broadening needs to be doubled in the effective operators for Z and X to achieve high accuracy between the evolution due to the effective master equation (S21) and the more complete master equation in Eq. (S5) (the exact prefactor in the W operators turns out to be less important). This can be attributed to interference between the blue-and red-detuned drives. Indeed, considering coherent excitation of a bright state consisting of both the blue-and the red-shifted dressed state suggests an increase of the broadening by a factor of two and thus the replacements of the effective detunings |∆ eff | 2 → |∆ eff | 2 +2nΩ 2 and |g eff | 2 → |g eff | 2 +2nΩ 2 . While these replacements are less rigorous, they are supported by our numerical simulations: It can be seen from Fig. 3 in the main manuscript that (1) the analytical scalings derived using the strong driving operators comprise upper bounds to the numerically obtained scalings and that (2) the effective operators for strong driving agree well with simulations of the more complete master equation.
Following the reasoning above, the power broadened decay rate for the Z pumping from n 1 to n 1 − 1 is found to be (S173) For the X depumping of |GHZ − we have and for the rate of the X process when acting on other states with n 1 (the "X toss" rate) Similarly, we find for the W scheme the "power-broadened" collective decay rate from |0 ⊗N to |W , (S176) For the individual decay process from |0 ⊗N to |W we have (S177) and for the individual decay process from |0 ⊗N to "as" (S178) Pumping of antisymmetric states with n 1 = 1 ("as") to |0 ⊗N , using only a single, resonant field, becomes As a consequence of the inclusion of terms for power broadening, increasing the driving strengths in the desired processes also increases the effect of power broadening in the desired decay rates. In the off-resonant decay rates γ eff ∝ γΩ g 2 the effect of power broadening is, on the other hand, negligible (given that Ω 2 g 2 ). The detrimental rates thus still increase for a growing Ω while the desired rates saturate. This is the limiting factor for the drive Ω which we will use to derive the possible Ω further down.

B. Population of the excited states
An additional reduction in the fidelity for strong driving comes from the population of the excited induced by the drive. In the following, we investigate this effect: The GHZ state, through its contribution from |1 ⊗N , is coupled to an excited state |ψ e = 1 √ N (|1..1e + |1..e1 + ... + |e..11 ). Despite being off-resonance, all tones of the driving field couple to the transition from |GHZ to the dressed states of |ψ e and |1 ⊗N |1 , with driving strengths of Ω For each tone the excited population is then approximately given by We also consider the excited population caused by X pumping. This requires representing the GHZ state in the X basis as in Eq.
(S109) and leads to an expression For W and A pumping we find the corresponding expressions These expressions are included in our numerical simulations of the effective dynamics in the strong driving regime and provide another limitation on Ω. Our analytical derivation of the scaling, on the other hand, involves upper bounds to certain error processes and turns out to take place in a parameter regime where the population of the excited states is not significant. Therefore, in the analytical considerations below we thus leave out its effect.

VII. SCALING ANALYSIS OF THE PREPARATION TIME AND THE ERROR OF THE GHZ PROTOCOL
In the following, we provide an analytical study of the scaling of the GHZ scheme. We derive an expression for the preparation time and optimize it by the choice of the parameters. The analysis is performed both for weak driving, beginning in Section VII A, and for strong driving, in Section VII E.

A. Optimization of the parameters for Z pumping alone
The different schemes presented in the main manuscript have in common, that the population of a nearly exponential number of states is pumped to subspaces with a polynomial number of states in a number of steps linear in the size of the system. For GHZ we have engineered strong decay processes from N − 1 ≥ n 1 ≥ 1 and for W from N ≥ n 1 ≥ 2. The decay is achieved using the 'Z' pumping which is only sensitive to n 1 and reduces this to n 1 − 1, finally leading to n 1 = 0 for GHZ and n 1 ≤ 1 for W state preparation. In order to assess the performance of the schemes it is thus important to know the time for a concatenated process consisting of many consecutive Z pumping steps. The rate of each individual decay is given by (see Eq. (S100)): The average time for this decay to occur is given by the inverse decay rate, For the total time required for pumping from an n 1 to an n 1 we add the average times for the intermediate steps We can also assign a total decay rate from an initial state with n 1 to a final state with n 1 , Γ n1→n 1 = 1 τ . It is however more useful to use the total preparation time and to minimize it by the choice of available parameters. Here, in particular the Rabi frequencies Ω (F ) Z of individual field tones F can be chosen, as well as the tunable decay rate γ e . Since the pumping occurs from a maximal n 1 of N − 1 to n 1 − 1, the worst case preparation time for Z pumping is (using Eq. (S186)) found to be given by The time from n 1 = 1 to |GHZ differs from that to n 1 = 0 by a factor of 2. This is due to the fact that only half of the population is pumped to |GHZ and the other half to |GHZ − , which is continuously depumped by the X pumping discussed below. Therefore, on average two attempts are required so that the preparation time is doubled, As anticipated already before Eq. (S107), for GHZ preparation we will from now on choose the parameter values κ b = κ c = 0, γ 0e = γ 1e = γ e /2, and abbreviate γ e ≡ γ. In order to obtain dimensionless optimization variables, we will furthermore write Ω (F =n1) Z =: A F Ω (for F = 1, 2, . . . , N − 1) with nonnegative dimensionless variables A F . The quantity Ω is a dimensionful frequency parameter, whose size has been chosen such that (for the weak driving calculation) all Ω ≤ Ω, i.e. such that the dimensionless parameters A F satisfy A F ∈ [0, 1]. In this notation, the above GHZ preparation time reads: where we have defined the function (S190) With the same abbreviations, the error rate from Z pumping alone reads, from Eq. (S107): where Our goal for now is to find parameters {A F } that minimize the Z-error for any given value of the GHZ preparation time (or, equivalently, minimize the GHZ preparation time for any fixed error value). Minimizing G({A F }) under the constraint H({A F }) ≡ H ≡ const by the method of Lagrange multipliers leads to Approximating the latter harmonic sum gives roughly η ≈ (log N )/ H, and we would have A F > 1 for some F (in particular, for F = 1) if η > 1/(N − 1), i.e. if H (N − 1) log N . However, one can still find the optimal assignment {A F } minimizing G({A F }) while obeying 0 ≤ A F ≤ 1 and H({A F }) ≤ H: where we defined the "ceil-1" function and η needs to be adjusted such that H({A F }) = H. The assignment (S194) means that A 2 F = 1 for F < N/(1 + 1/η) and A 2 F = η(N − F )/F for F ≥ N/(1 + 1/η). That (S194) is the unique optimal solution can be checked by the Karush-Kuhn-Tucker (KKT) conditions [52, Section 5.3.3], using that the function G({A F }) to be minimized and the constraint function H({A F }) are both strictly convex in their arguments. For η ∈ [1/(N − 1), N − 1] (such that the selection between the two cases in (S195) happens at some F ∈ [1, N − 1]) one can thus compute: and For η = const > 0, the error in both estimates is O(1) as N → ∞ and thus irrelevant compared to the log N terms. In particular, if one wants the numerical factor H in the preparation time (S189) to scale like log N (like optical pumping), then one needs η ≥ O(1), which we achieve by letting η = const as N → ∞.
The error E can generally be obtained by comparing the preparation time of the desired state of the protocol and the loss rate out of it. Using Eqs. (S189) and (S191), the error E is found to be given by Thus, to achieve a desired error E (which may be N -dependent, i.e. E = E(N )), one can adjust γ appropriately. If we assume a relation Ω = αγ in order to limit Ω to the weak driving regime (and where the number α may or may not depend on N , i.e. α = α(N )), we can plug this back into the worst-case preparation time (S189) to obtain where the last estimate holds for η = const as N → ∞ and we have neglected lower-order terms in N . The prefactor is minimized for η = 2, leading to a minimal preparation time (for the desired error E): To summarize the optimal parameter choices for the scenario considered here: 5. Four-compartment-model of the process creating |GHZ .
• The A F for 1 ≤ F ≤ N − 1 have to be chosen as follows (cf. Eq. (S194) with the optimal choice η = 2): i.e. A F = 1 for 1 ≤ F ≤ 2N/3, and A F = 2(N − F )/F for 2N/3 < F ≤ N − 1. This leads to: • The choice of γ ≡ γ e , and consequently of Ω ≡ αγ and γ 0e = γ 1e = γ/2, is given by: We also set κ b = κ c = 0. This leads to: The appearance of N 1/2 in the scaling of τ n1=N −1→GHZ can be traced back to the fact that the noise Γ − acts on each of the N atoms (i.e. the prefactor of N in Eq. (S107)).

B. Compartment model and effective rates
While so far we have only discussed the Z pumping, we will now make a simplified model of both processes, Z and X, that create the |GHZ state. This model is more pessimistic w.r.t. the preparation time and the error treatment than the actual (Lindbladian) dynamics described in Section IV. Nevertheless, this model will still yield a scaling of the preparation time as τ GHZ ∼ (N 1/2 log 2 N )/ √ E with the number of qubits N and the error E, just like Eq. (S207), which thus shows that this is indeed the best achievable scaling. Our simplified model is shown in Fig. 5 and consists of four compartments as we explain now. We split the 2 N -dimensional Hilbert space into the subspaces with n 1 = 1, 2, . . . , N − 1, counting the number of |1 -states appearing in a computational basis state (as in Section IV). All states with n 1 = N − 1 are put into compartment 1, as they are furthest away (in pumping time) from the GHZ states; the states with n 1 = 1, . . . , N − 2 are then put into compartment 2. The two Hilbert space dimensions belonging to n 1 = 0 and to n 1 = N are split into the compartments 3 and 4 of Fig. 5, corresponding to the |GHZ − and the |GHZ states, respectively. The desired Z pumping process creates each of those GHZ states with equal rate 1 2 Γ + Z out of compartment 2 where using Eq. (S186) we find Here, we have again plugged in γ 0e = γ/2, κ b = 0 and the optimal A j from Eq. (S208) and followed the same computation as for (S199), neglecting lower-order terms. Similarly, the rate from compartment 1 to compartment 2 is given by Eq. (S184): As compartment 1 is the furthest away from the desired GHZ state, we pessimistically model the X-toss (i.e. the action of the X process on the states other than |GHZ and |GHZ − , see Eq. (S131)) to throw any state back to the n 1 = N − 1 compartment, with a rate Γ toss X to be computed below. Similarly, we model the errors Γ − Z and Γ − X affecting |GHZ such that they move the |GHZ -state back to compartment 1. By the same rationale, even the "good" X process Γ + X is modelled to take the un-wanted state |GHZ − back to compartment 1. The detrimental Z rate for the compartment model as computed with the parameters from the previous subsection, Eq. (S208) (see also Eqs. (S191), (S192), and (S203) with the optimal η = 2) is For the X rates, it was found in Section IV B that only the Rabi oscillations Ω (F ) X with odd index F should be turned on. Similar to Section VII A, we write Ω X Ω with the dimensionless parameters A (F ) X which we discuss below. The "good" X rate Γ + X for the compartment model of Fig. 5 is then given by Eq. (S127): Since the (normalized) probability distribution { N F /2 N −1 } N F =1,3,... is strongly peaked around the values F ∼ N/2 (similar to the full binomial distribution), the exact functional choice of the X coefficients A (F ) X (for odd F ) does not really matter in the limit of large N , as long as it is not exponentially fine tuned. Since a similar binomial distribution occurs in the X error rates below as well and as the same reasoning applies, we can take all A (F ) X to be equal as a very good approximation: With this, the rate Γ + X evaluates to (with exponentially good accuracy for large N ): As we will see below, with this choice the error induced by the X error rate is smaller than the one from the Z process and thus not a major limitation to the preparation of the entangled state. The X error rate Γ − X is given by (S130), and we simplify it again with our above parameter choices, making approximations of the binomial distribution and the other sum which are good in the large-N limit: In the expression above, [N/2] even denotes the next higher integer number of N/2. Having the limit of large N in mind, we evaluate the last sum as follows: neglecting subleading terms in N . This finally gives: Finally, the X toss rate Γ toss X in Fig. 5 is given by Eq. (S131), which with our parameter choices becomes [note that (S131) is half of the "good" rate (S127), which we have computed in Eq. (S219) already]: From the process in Fig. 5 one can see that the optimal parameters A X have to be chosen such that the good X-process and the good Z-process have about the same rate: If the X pumping rate Γ + X is too weak, population will accumulate in |GHZ − by the Z pumping (Γ + Z ). On the other hand, a too strong X pumping will hinder the preparation mechanism through the X toss effect. We thus set the rates for the desired processes, Z pumping and X depumping to be equal, This results in: (note that, for all N ≥ 2, this choice is consistent with the requirement A X ≤ 1 for the weak-driving analysis). With these choices made, we summarize the parameters and effective rates for the four-compartment model of Fig. 5 found so far: • For the Z-pumping, we make the choice of parameters found to be optimal in Eq. (S208): meaning that Ω • For the X-pumping we take (see Eq. (S231)): • Then one obtains for the effective rates in Fig. 5: Γ toss • The total error rate (leading from compartment 4 to compartment 1) is thus 9N log 2 N 16 1 + 5 9 log 2 N . (S239)

C. Transition matrix, stationary error, and GHZ preparation time
The transition matrix for the 4-compartment model of Fig. 5 is: The steady-state population p ∞ := (P 1 (∞), P 2 (∞), P 3 (∞), P 4 (∞)) is given as the solution (normalized to the sum of entries being 1) of the equation T p ∞ = 0. This gives: The steady-state fidelity is just F = P 4 (∞), and the error is thus (Here the approximation in the second line was made for analytical convenience and gives a slightly pessimistic bound.) In the scaling with large N , this expression for E agrees with the one implied by Eq. (S211) that was found by other means before, and the prefactor is similar. Thus, to achieve a desired stationary error E, we need to adjust γ such that: Below we will use this expression instead of the results obtained in (S211) which were derived by considered only the Z pumping. So far we have discussed the Z pumping separately from the X pumping, deriving an individual characteristic time τ n1=N −1→GHZ . It now remains to derive an analytical expression for the total GHZ pumping time that is obtained in the presence of X pumping, using the parameters (S244) and (S208). In order to factor out the dependence of the stationary error E and to obtain a tractable analytical expression, we make (for the computation of the GHZ preparation time) an approximation to the transition matrix (S240) by dropping the small terms leading out of the GHZ state (these terms vanish in the limit of E → 0). That is, we set the fourth column in (S240) to zero: For the initial population vector (P 1 (0), P 2 (0), P 3 (0), P 4 (0)) = (0, 0, 1, 0), which corresponds to the whole population being in the worst state |GHZ − of Fig. 5, we have the following evolution: Thus, on time-scales larger than any fixed t 0 , the transition from the worst state |GHZ − to the desired state |GHZ happens at least as fast as in an exponential decay process, with an approximation Γ + for the effective exponential rate computed as: Note that the fraction in the last expression will depend both on (Γ + Z t 0 ) and on N , since the N -dependence of the transition matrix T + in Eq. (S245) cannot be factored out completely. To get a meaningful expression for the exponential rate, the timescale t 0 should be chosen comparable to the other relevant timescales of the process. For definiteness we will thus set t 0 := 1/Γ + Z throughout. The characteristic preparation time of |GHZ in the sense of an exponential rate is then: In the limit of large N , the square root inside the bracketed expression will tend to 1. And also the fractional expression involving P 4 (t 0 ) will tend to a constant number independent of N , since large values of the first column in the matrix in (S245) mean that the transition time out of the first compartment is insignificant compared to the other transition times, which are all independent of N . This is also easily seen numerically. In the following table, we evaluate the factor in square brackets (which we call b(N )) in Eq. (S250) for different values of N : The GHZ preparation time is then: where b(N ) tends to about 13 as N → ∞ (see table above). The parameter choices for this can be found in Eqs. (S232), (S233), and (S244) together with Ω = αγ, γ 0e = γ 1e = γ, κ b = κ c = 0.

D. Analysis of the "dynamical problem"
Here we analyze the "dynamical problem". This means that, for a certain exponential form of the time-evolution of the GHZ error E(t) (or, equivalently, of the GHZ fidelity F (t) = 1 − E(t)) in an effective model, we compute and minimize the time t it takes to achieve a desired target error E. The main result will be that, up to an additional factor of log(1/E), the scaling of the preparation time with the particle number N and with the error E is the same as in the previous analytical approaches, see e.g. Eqs. (S212) and (S251) (with E replaced by E).
In the following, we write again Ω = αγ to be able to limit our treatment to the weak-driving regime by suitable choices of the number α. Furthermore, to keep the main derivation as general as possible, we write the rates into and out of the desired GHZ state as with functions f ≡ f (N ) and h ≡ h(N ). Later, we will evaluate our results for the functions which is motivated by Eq. (S239), and which is chosen such that the ratio Γ − /Γ yields the stationary error from Eq. (S243). Finally, we make the following basic ansatz for the time evolution of the error: Note that in the limit of t → ∞ the error indeed converges to Γ − /Γ with an exponential rate given by κΓ, where κ > 0 can be a dimensionless constant to adjust the effective decay rate in a model where the effective decay rate κΓ does not match with the steady state error Γ − /Γ (see Eqn. (S278) and below). Since we are interested in the regime of small stationary error Γ − /Γ, we can approximate and continue: where we abbreviate with a constant c and a "rescaled time" τ as follows: We can treat c as a free optimization variable, since γ is a freely adjustable parameter, and thus we can adjust c to any nonnegative real number by choosing γ appropriately (even if g and h(N ) and f (N ) are fixed). Furthermore, τ is essentially the same as the physical time t, but rescaled by a fixed number (which in particular depends on N and g).
The optimal time τ can be obtained by plugging this expression into (S262) and simplifying the expression, but there is a less direct and somewhat easier way: Observe from Eqs. (S262) and (S263) that where in the last step we used the asymptotic approximation of our branch of the Lambert W function: W −1 (x) = log(−x) − log(− log(−x)) + O(1) as x → 0. This is justified when we are interested in the case of very small or asymptotically vanishing error E → 0.
Finally, we can use the above expressions to solve Eqs. (S260) and (S261) for the physically interesting optimal parameters γ = γ(N, E) and t GHZ = t(N, E) using the values given in Eqs. (S268) and (S270). When we use the choices for f (N ) and h(N ) given in Eqs. (S254) and (S255), we obtain: Note that for small dynamical error E we have W −1 → 0, such that Eq. (S273) approaches the previous result in Eq. (S244). This is due to the fact that for small E the stationary part of the error dominates. The GHZ preparation time (defined as the time to reach a GHZ error of value E) is: The approximations of the Lambert W function used above, i.e. 1 + 2/W −1 → 1 and W 2 −1 + 2W −1 → log(1/E), are good only for quite small E. For E = 0.1, one should instead use the exact value W −1 (−2E/e 2 ) = −5.27, leading to 1 + 2/W −1 = 0.788 and W 2 −1 + 2W −1 = 4.15, which makes that Eq. (S275) is by a factor 1.8 lower than Eq. (S274). For E = 0.03 the corresponding value is W −1 (−2E/e 2 ) = −6.72, leading to 1 + 2/W −1 = 0.838 and W 2 −1 + 2W −1 = 5.63, which makes that Eq. (S275) is by a factor 1.6 lower than Eq. (S274). We can also find the relation between the stationary error E = Γ − /Γ (from Eq. (S256)) and the dynamical error E: which is independent of N . Thus, using the above values, for E = 0.1 we get for the stationary error at the optimal parameters E = 0.62E, whereas for E = 0.03 we get E = 0.70E. It now remains to fix the value of κ (appearing in (S274)-(S275)) in an appropriate way, namely such that the model (S256) matches the GHZ preparation process from Sections VII B and VII C as well as possible. For this, note that the effective GHZ preparation rate Γ GHZ from Section VII C can be inferred from Eq. (S248): Using the value of Γ from Eqs. (S252) and (S255), we can solve for κ: Evaluating this as in Section VII C, we get the following table: When taking power broadening into account (see Section VI A), then instead of Eq. (S184) from the weak driving scenario, the favorable transition rates of the Z process are now given by Eq. (S173) (for n 1 = 1, 2, . . . , N − 1): where we have again used the parameter values and abbreviations γ e ≡ γ, γ 0e = γ/2, Ω (F =n1) Z = Ω F as in Section VII A. Thus, instead of (S188), the Z pumping time is now: The error rate from (S191) remains unchanged: As below (S192), we now minimize Γ GHZ→?,γ,Z (as a function of the variables {Ω f }) while keeping τ n1=N −1→GHZ constant. This leads, by the method of Lagrange multipliers, to: where λ is a Lagrange multiplier (that has units of frequency 2 ), which we will determine later. Plugging this back into (S282) and (S283), we get: Thus, the stationary Z error is (cf. (S204)): Thus, when the desired stationary Z error E Z is given, then λ and γ are determined by each other (since g and N are fixed): Plugging this back into (S285), we get: For fixed N , g, E Z , and γ, this is the minimal Z pumping time (i.e. minimized over all choices of pumping rates {Ω f }). Since we do not want to fix γ a priori, we minimize the last expression over γ, finding that the optimal choice is The corresponding pumping strengths can now be computed from (S284) and (S288): and the optimal Z pumping time for given N , g, E Z is then (compare to (S207) without power broadening): From now on we set (N − 1) N , as this neglects only subleading terms.
We also take power broadening into account for the desired X pumping rate (see (S174)) and for the X toss rate (see S175), which we both compute as in (S216)-(S219) and (S229), and we take the detrimental X rate from (S228): Γ toss Now we adjust the X parameters Ω and γ f such that the X rates agree with the corresponding Z rates (cf. Section VII B), i.e.
gE 3/2 Z N 3/2 log N . To solve this for γ f and Ω exactly, one would have to solve a cubic equation. As this is quite cumbersome and uninformative, we are looking for solutions which satisfy the following scaling ansaetze: γ f N α (log N ) β and Ω N φ (log N ) ψ . Then, one finds actually two possible solutions for Ω and γ f , which lead to the desired scaling of Γ + X and Γ − X (as N becomes large). One of the solutions is given by Formally, another solution exists, but this is in the extremely saturated regime where our effective operators do not apply.
Finally, we need to find the relation between the above rate Γ + Z := 1/(τ N −1→GHZ,Z ) and the total GHZ preparation time τ GHZ ≡ 1/Γ + (which includes the errors and X toss) and the total stationary error E, just as in Section VII C. For this, we consider a simplified 3-compartment model which constitutes a very good approximation in the large-N -limit, in which the parameters (S297) were computed in the first place. Here, compartment A comprises the states with n 1 = 1, 2, . . . , N − 1, compartment B is the state |GHZ − and compartment C the state |GHZ . Then the transition matrix is (cf. (S240) and Fig. 5): again with Γ − = Γ − X + Γ − Z , and again as in Section VII C we denote by T + the transition matrix without the "bad" rates Γ − . From the stationary vector p ∞ satisfying T p ∞ = 0 we can again compute the total stationary error E: We therefore have to set E Z = E/5 in all previous expressions. The relation between Γ + Z and Γ + ≡ 1/τ GHZ is similar to Eq. (S247): where again P C denotes the population in the GHZ state when starting from |GHZ − state. Computing this for the above transition matrix and plugging in the values for Γ + Z and E from above, we obtain: Thus, in the large-N limit, the final GHZ preparation time τ GHZ is: This is achieved by the following parameter choices in terms of g, N , E: Again, as in Eq. (S277), if one is interested in the dynamical error E instead of the static error E, one should everywhere set E = 0.62E for E = 0.1 (or E = 0.70E for E = 0.03). Furthermore, the GHZ preparation time is then prolonged by an additional factor of log(1/E) (see Eq. (S275)).  6. 5-compartment-model of the W process. The "good" transition rates are shown as bold lines with arrows, whereas the error processes are shown with dashed lines.

VIII. SCALING ANALYSIS OF THE PREPARATION TIME AND THE ERROR OF THE W PROTOCOL
Similar to our scaling analysis of the GHZ preparation time in Section VII, we analyze the scaling of the W preparation time. As already anticipated from the transition rates computed in Sections V and VIII B, we will again employ a transition rate model, now consisting of five compartments as shown in Fig. 6.

A. Weak driving analysis
Compartment 5 encompasses all states with n 1 ≥ 3 excitations. As can be seen from Section V, none of the error processes leads into this compartment, such that this compartment is continually emptied. As can be seen from Fig. 6, the time to reach the W state (compartment 1) is longest time when starting from n 1 = N , which we will henceforth assume as the worst case. Thus, the transition rate Γ + 54 from compartment 5 into compartment 4 (consisting of the states with n 1 = 2 excitations) is computed from the N − 2 individual transitions N → (N − 1) → . . . → 3 → 2. With Eq. (S139) and following the rule of inverse addition of decay rates we get: where we have made the parameter choice that all Ω (F =n1) Z for n 1 = 3, . . . , N are equal to each other and we call their value Ω (F ≥3) Z . The inverse of the sum on the right-hand-side equals (log N ) −1 , up to a correction that is smaller than (log N ) −2 , i.e. of lower order. As physically and implementationally reasonable parameter choices (see also Section V), we will throughout set Furthermore, we introduce a dimensionless constant β to write κ b as: Later on, it will be seen that β is to be chosen as a constant (of order 1) in order to obtain an optimal scaling of the W preparation time. From Eq. (S139) it can be seen that the desirable decay rate saturates for too strong drives since the transition is saturated. Hence, driving stronger than Ω (F ≥3) Z ≈ γ e + κ b = γ e (1 + β) does not significantly increase the preparation rate, but only the loss rates. We thus re-parametrize Ω (F ≥3) Z with a dimensionless constant α n : With these substitutions and assumptions, the transition rate Γ + 54 from (S308) reads: With the same substitutions, Eq. (S138) yields the transition rate from compartment 4 to compartment 3 (consisting of the states with n 1 = 1 excitation and outside of the sector with permutation symmetry, i.e. the "antisymmetric" n 1 = 1 states): where the dimensionless constant α Z (potentially different from α n ) has been introduced: Here we have made the worst case assumption that decay from states with n 1 = 2 leads to the antisymmetric states in n 1 = 1. The favorable transition rate from compartment 3 to compartment 2 (consisting of the n 1 = 0 state) is given by Eq. (S173). As the parameter Ω A does not appear in the loss rates (Eqs. (S153)-(S155)) in the present approximation, it will clearly be best to choose it as large as possible. However, when we include power broadening, the decay rate saturates at around Finally, the transition rate from compartment 2 to compartment 1 (the desired W state) is given by Eq. (S165). Again due to the restrictions set by power broadening, we must have Ω W ≤ γ e + κ b = (1 + β)γ e . We reparametrize with a dimensionless parameter α W : (Note, here we take the factor √ N into the parametrization for notational convenience, since in both the good rates and the loss rates the parameter Ω W always appears in the combination N Ω 2 W , see Sections V and VIII B). Thus, we get from Eq. (S165): For the full W preparation time τ W , we get from Eqs. (S313), (S314), (S316), (S318) and the rule of the inverse addition of decay rates: Next we compute the loss rates from the W state, i.e. the rates corresponding to the dashed lines in Fig. 6. For this, observe that the parameter κ c occurs only in the loss rates (Section V B), but not in the gain rates; thus, it will be optimal to set κ c = 0 .
Given the gain and loss rates just computed, the intensity matrix (Markov transition kernel) for the compartment model of Fig.  6 is: between the decay rates that is independent of N and has a value r ∼ 1. With this substitution, the stationary error becomes: (1 + β) 2 (1 + 2β) 4 .

(S335)
Note that the expression in the square brackets is dimensionless, as are all parameters occuring in it. We also assume that g is fixed from the beginning. Thus, for any given scaling of the dimensionless parameters (which we will determine below), we must choose a small enough γ e in order to achieve the stationary error E. To get E, we must choose where [brackets] denotes the expression in square brackets in Eq. (S335). Plugging this into the W preparation time (S320) (using (S334), so that all quantities except γ e are dimensionless), we obtain: Here, there are terms log(N ) which grow with the qubit number, while the terms with r are constant and can thus be omitted for large N . Given that α Z and α W enter in the same way, we can assume that there is little gain in having them different. We thus set them equal to each other. We then obtain for the preparation time Eq. (S338) identifies the resulting scaling in the weak driving regime. Similar to the result for GHZ it goes as the inverse of the strength of the driving, parametrized by α. To investigate the true scaling we thus have to perform a strong driving analysis. This is done below, but in short we find that the main limitation comes from α n . This is due to the fact that α n enters in Γ + 54 which describes the rate at which we can depump the manifold with n 1 ≥ 2. It necessarily involves O(N ) decays and a resulting time O N γ . Comparing to Eq. (S313) we see that this puts a restriction α n ∼ 1 √ N . Inserting this in Eq. (S338) gives us the scaling which only differs from the result obtained below by logarithmic corrections.

B. W scaling analysis for strong driving
We will analyze here the scaling of the W preparation time using effective operators accounting for the power broadening effect (Section VI A). As in Section VIII A, we employ again the transition rate model from Fig. 6. To keep the analysis manageable, we will again assume simple and physically well-motivated relations between the decay rates and express all of them in terms of a single parameter, which we call γ throughout (cf. also Section VIII): γ 0e = γ 1e = 1 2 γ e = 1 2 γ, κ c = 0.
From Eq. (S173), we get for the rates of the Z pumping (for n 1 = 2, 3, . . . , N ): where we denote the different pumping strengths for the optical Z pumping by Ω F with F = 2, 3, . . . , N . By the law of inverse addition of decay rates, we thus obtain the following "good" rates Γ + 54 and Γ + 43 (cf. Fig. 6): The only errors caused by the Z pumping drives Ω F are, by Eqs. (S140) and (S141): Besides in these errors, the Z pumping drives Ω F appear in the preparation time τ W = (Γ + 54 ) −1 + (Γ + 43 ) −1 + (Γ + 32 ) −1 + (Γ + 21 ) −1 (see below and also e.g. Eq. (S319)), and they appear only in the combination (Γ + 54 ) −1 + (Γ + 43 ) −1 = 2(N −1) To find the optimal choice of {Ω F } N F =2 , we therefore want to minimize the total Ω F -error in Eq. (S346) for each constant value of (Γ + 54 ) −1 + (Γ + 43 ) −1 , or equivalently, for each constant value of . The method of Lagrange multipliers thus leads us to consider the following Lagrangian: where λ is a Lagrange multiplier that has dimensions of [frequency] 2 . Solving the minimizer conditions dL dΩ F = 0 for all F = 2, . . . , N , we get the solutions Thus, the optimal choice for the (N − 1) parameters {Ω F } depends on only one parameter λ, which we will choose below in the optimal way. Plugging (S348) back into (S345), into (Γ + 54 ) −1 + (Γ + 43 ) −1 , and into the error rates (S346), we get: Γ W→0,γ0,Z = Γ W→as,γ1,Z = γ 4g 2 The good rate Γ + 32 is, from Eq. (S177): Due to the choice κ c = 0, the parameter Ω A appears only in this good rate, but not in any error rate (see Section V B). Thus, Ω A should be chosen in such a way as to maximize Γ + 32 , which happens for Ω A γ f . As we have assumed g to be the largest parameter in our derivation, we set, however, Ω A = γ f . We thus have: By Eq. (S176), the rate Γ + 21 (with collective decay) is: We examine now the exact optimal N -scaling in the asymptotic regime. Note that the expression under the square root in (S371) satisfies: Thus, the preparation time τ W in (S371) satisfies: And indeed, such a N -scaling like N (log N ) 1/4 can be achieved in (S371) by choosing α = α 0 and β = β 0 /(log N ) 1/2 with positive constants α 0 and β 0 . Then the overall asymptotic prefactor is minimized for α 0 = β 0 / √ 6 and small positive constant β 0 . Note that, due to the log-term in Eq. (S371), the limit β 0 → 0 is not, however, sensible for any finite N . For this reason, and since the numerical prefactor can never become smaller than 8 √ 6 = 4.42 anyway, we may for example want to choose β = 1/(log N ) 1/2 and α = 1/ √ 6, in which case the asymptotic W preparation time is: For the simpler choice α = β = 1, the scaling becomes only a little worse: Note that, due to the log-terms, the asymptotic N -scaling of the preparation time (S371) is relevant only for rather large N (e.g. only when 2 √ 6(log N ) 1/2 > (4 + 17α)), and one should rather use the parameter values provided in the preceding table. For realistic (small) values of N , these parameter choices lead to somewhat better τ W than any of the choices that led to Eqs. (S376) and (S377). To conclude, we summarize the above parameter choices: • choose the dimensionless parameters α, β such that, for the given N , they minimize the expression (S371), as done in the above table for N = 2, . . . , 10. (Alternatively, choose for example α = β = 1).
Similar as in the GHZ scheme, the dynamical error E is obtained by multiplying the static error E with a factor (cf. Eq. (S277), e.g. E = 0.62E for E = 0.1 or E = 0.70E for E = 0.03. The W preparation time is then prolonged by an additional factor of log(1/E) (see Eq. (S275)).