Autonomous multipartite entanglement engines

TL;DR

This paper introduces autonomous entanglement engines that generate various multipartite entangled states using only thermal baths and static interactions, advancing practical quantum information processing.

Contribution

It proposes a general architecture for heralded generation of multipartite entangled states using incoherent thermal resources, with conditions for specific states like GHZ, Dicke, and cluster states.

Findings

Able to generate GHZ, Dicke, and cluster states

Uses only thermal baths and static interactions

Provides conditions for target state generation

Abstract

The generation of genuine multipartite entangled states is challenging in practice. Here we explore a new route to this task, via autonomous entanglement engines which use only incoherent coupling to thermal baths and time-independent interactions. We present a general machine architecture, which allows for the generation of a broad range of multipartite entangled states in a heralded manner. Specifically, given a target multiple-qubit state, we give a sufficient condition ensuring that it can be generated by our machine. We discuss the cases of Greenberger-Horne-Zeilinger, Dicke and cluster states in detail. These results demonstrate the potential of purely thermal resources for creating multipartite entangled states useful for quantum information processing.

Equations

$$H_{\text{free}}=\sum _{k=1}^N{H}_k=\sum _{k=1}^N\left(\sum _{l=1}^2{\Delta }_k^{(l)}|l{⟩}_k⟨l|\right).$$

$$\frac{d\rho }{dt}=-i[H_{\text{free}}+H_{\text{int}},\rho ]+\mathcal{L}(\rho ).$$

$$\mathcal{L}(\rho )={\mathcal{L}}_k(\rho )=\sum _{k=1}^N{\gamma }_k({\tau }_k{\otimes }_k{\mathrm{Tr}}_k(\rho )-\rho ).$$

$${\rho }^{\prime }=\frac{\Pi {\rho }_{\infty }\Pi }{\mathrm{Tr}\left({\rho }_{\infty }\Pi \right)}$$

$$H_{\text{int}}=g\left(|\bar{\psi }⟩⟨R|+|R⟩⟨\bar{\psi }|\right),$$

$$|\psi ⟩=\sum _{\mathbf{n}\in S_{\psi }}{c}_{\mathbf{n}}|\mathbf{n}⟩,$$

$$\displaystyle\frac{1}{2}\sum_{k=1}^{N}\bigg{[}R_{k}n_{k}\Delta^{(1)}_{k}$$

$$\displaystyle-\frac{1}{2}\sum_{k=1}^{N}\bigg{[}R_{k}\Delta^{(2)}_{k}\bigg{]}=0,$$

$$n_{k^{\prime}}\Delta^{(1)}_{k^{\prime}}+\sum_{k\neq k^{\prime}}\bigg{[}(1-n_{k})\Delta^{(1)}_{k}+n_{k}\Delta^{(2)}_{k}\bigg{]}-\Delta^{(2)}_{k^{\prime}}=0.$$

$$H_h$$

$$H_c$$

$$H_{\text{int}}=g(|2\bar{0}⟩⟨1\bar{1}|+|2\bar{0}⟩⟨0\bar{2}|+|1\bar{1}⟩⟨2\bar{0}|+|0\bar{2}⟩⟨2\bar{0}|),$$

$$|D_l^N⟩=\frac{1}{\sqrt{\left(\genfrac{}{}{0ex}{}{N}{l}\right)}}\sum _s{\sigma }_s\left[|1{⟩}^l\otimes |0{⟩}^{N-l}\right],$$

$${\Delta }_h^{(1)}$$

$${\Delta }_h^{(2)}$$

$$|C\rangle=\frac{1}{2}\big{(}|0110\rangle+|0101\rangle+|1010\rangle-|1001\rangle\big{)}.$$

$$H_h$$

$$H_c$$

$$|\psi ⟩=\sum _{\mathbf{n}\in S_{\psi }}{c}_{\mathbf{n}}|\mathbf{n}⟩$$

$$P_o\left(2\frac{{\gamma }_h}{3}+{\gamma }_c(N-1)\right)=\frac{{\gamma }_h}{3}\left(P_{\bar{\mathbf{n}}}+P_{o^{\prime }}\right),$$

$$\frac{P_o}{P_{\bar{\mathbf{n}}}}=\frac{{\gamma }_h}{3(N-1){\gamma }_c+{\gamma }_h}.$$

$$\sum _{o\in S_{\bar{\psi }}}{\bar{P}}_o+\sum _{o\in G_h}{\bar{P}}_o=1.$$

$${\bar{P}}_{\text{S}}\left(1+\frac{{\bar{P}}_{\text{G}}}{{\bar{P}}_{\text{S}}}\right)=1$$

$${\bar{P}}_{\text{S}}=\left(1+\frac{{\bar{P}}_{\text{G}}}{{\bar{P}}_{\text{S}}}\right)={\left(1+\frac{{\gamma }_h}{3(N-1){\gamma }_c+{\gamma }_h}\right)}^{-1}.$$

$${\bar{P}}_{\bar{\mathbf{n}}}=⟨\bar{\mathbf{n}}|{\rho }^{\prime }|\bar{\mathbf{n}}⟩=|c_{\bar{\mathbf{n}}}{|}^2.$$

$$0={\dot{\rho }}_{n,n^{\prime }}=-i⟨n|[H,\rho ]|n^{\prime }⟩+\frac{{\gamma }_h}{3}⟨n|\text{openone}\otimes {\mathrm{Tr}}_1\left(\rho \right)|n^{\prime }⟩+\sum _{k=2}^N{\gamma }_c⟨n|\left(|0⟩⟨0|{\otimes }_k{\mathrm{Tr}}_k\left(\rho \right)\right)|n^{\prime }⟩-\left({\gamma }_h+{\gamma }_c\right){\rho }_{n,n^{\prime }}.$$

$$\begin{gathered} ⟨n|[H,\rho ]|n^{\prime }⟩=g⟨n|\left(|\bar{\psi }⟩⟨R|+|R⟩⟨\bar{\psi }|\right)\rho -\rho \left(|\bar{\psi }⟩⟨R|+|R⟩⟨\bar{\psi }|\right)|n^{\prime }⟩ \\ =g\left(⟨n|\bar{\psi }⟩⟨R|\rho |n^{\prime }⟩+⟨n|R⟩⟨\bar{\psi }|\rho |n^{\prime }⟩-⟨n|\rho |\bar{\psi }⟩⟨R|n^{\prime }⟩-⟨n|\rho |R⟩⟨\bar{\psi }|n^{\prime }⟩\right). \end{gathered}$$

$$⟨n|[H,\rho ]|n^{\prime }⟩=\{\begin{array}{ll}g\left(c_n{\rho }_{R,n^{\prime }}-c_{n^{\prime }}^{\ast }{\rho }_{n,R}\right)& \text{if }n,n^{\prime }\in S_{\bar{\psi }}\\ 0& \text{if }n,n^{\prime }\notin S_{\bar{\psi }}\text{ and }n,n^{\prime }\ne R\end{array}.$$

$$⟨n|\text{openone}\otimes {\mathrm{Tr}}_1\left(\rho \right)|n^{\prime }⟩={\delta }_{n_1,n_1^{\prime }}\sum _j{\rho }_{j\bar{n},j{\bar{n}}^{\prime }},$$

$$⟨n|\left(|0⟩⟨0|{\otimes }_k{\mathrm{Tr}}_k\left(\rho \right)\right)|n^{\prime }⟩={\delta }_{n_k,0}{\delta }_{n_k^{\prime },0}\sum _{j_k}{\rho }_{\overleftarrow{n}{j}_k\overrightarrow{n},{\overleftarrow{n}}^{\prime }{j}_k{\overrightarrow{n}}^{\prime }},$$

$$0={\dot{\rho }}_{n,n^{\prime }}=-ig\left(c_n{\rho }_{R,n^{\prime }}-c_{n^{\prime }}^{\ast }{\rho }_{n,R}\right)+\frac{{\gamma }_h}{3}{\delta }_{n_1,n_1^{\prime }}\sum _j{\rho }_{j\bar{n},j{\bar{n}}^{\prime }}-\left({\gamma }_h+{\gamma }_c\right){\rho }_{n,n^{\prime }}.$$

$$0={\dot{\rho }}_{n,n^{\prime }}=-ig{c}_n{c}_{n^{\prime }}^{\ast }\left(L^{\ast }-L\right)+\frac{{\gamma }_h}{3}{\delta }_{n_1,n_1^{\prime }}{\rho }_{n,n^{\prime }}-\left({\gamma }_h+{\gamma }_c\right){\rho }_{n,n^{\prime }}$$

Full text

*Introduction.—*Quantum thermal machines combine quantum systems with thermal reservoirs at different temperatures and exploit the resulting heat flows to perform useful tasks. These can be work extraction or cooling, in analogy with classical heat engines and refrigerators, but may also be of a genuinely quantum nature. In particular, it is possible to devise entanglement engines – thermal machines generating entangled quantum states. Entanglement is a key resource for quantum information processing but is generally very fragile and easily destroyed by environmental noise. It is nevertheless possible to exploit dissipation to create and stabilise entanglement Plenio et al. (1999); Plenio and Huelga (2002); Schneider and Milburn (2002); Kim et al. (2002); Jakóbczyk (2002); Braun (2002); Benatti et al. (2003); Hartmann et al. (2006); Quiroga et al. (2007); Burgarth and Giovannetti (2007); Kraus et al. (2008); Diehl et al. (2008); Verstraete et al. (2009). This was studied in a variety of settings and physical systems Cai et al. (2010); Kastoryano et al. (2011); Žnidarič (2012); Bellomo and Antezza (2013); Reiter et al. (2013); Schuetz et al. (2013); Walter et al. (2013); Ticozzi and Viola (2014); Boyanovsky and Jasnow (2017); Hewgill et al. (2018); C. K. Lee (2019) and dissipative entanglement generation using continuous driving was experimentally demonstrated, mainly for bipartite states Krauter et al. (2011); Barreiro et al. (2011); Shankar et al. (2013); Lin et al. (2013).

Autonomous entanglement engines represent a particularly simple case. Here, entanglement can be generated dissipatively with minimal resources, using only time-independent interactions and contact to thermal reservoirs at different temperatures. No driving, coherent control, or work input is required. For the bipartite case, a two-qubit entangled state can be generated in a steady-state, out-of-thermal-equilibrium regime Brask et al. (2015). Although the entanglement produced by such machines is typically weak, it can be boosted via entanglement distillation Bennett et al. (1996), or by coupling to negative-temperature Tacchino et al. (2018) or joint baths Man et al. (2019). In fact, applying a local filtering operation to the steady state of a bipartite entanglement engine can herald maximal entanglement between two systems of arbitrary dimension Tavakoli et al. (2018a).

These first results show that using dissipative, out-of-equilibrium thermal resources offers an interesting perspective on entanglement generation. A natural question is whether this setting could also be used to generate more complex forms of entanglement, in particular entanglement between a large number of subsystems. It is of fundamental interest to understand the possiblities and limits of thermal entanglement generation. In addition, such multipartite entangled states represent key resources, e.g. for measurement-based quantum computation, quantum communications, and quantum-enhanced sensing and metrology. The creation and manipulation of complex entangled states is therefore of strong interest for many experimental platforms, although typically very challenging in practice.

Here, we propose autonomous entanglement engines as a new route to the generation of multipartite entanglement and explore their potential. A first question is, which types of multipartite entangled states can be created. We present a sufficient condition for a given target $N$-qubit state to be obtainable. Specifically, for any target state satisfying our criterion, we construct an autonomous entanglement engine that will generate this state. The engine consists of $N$ interacting qutrits (three-level systems), each qutrit being locally connected to a thermal bath. From the resulting steady state, a local filtering operation then leads to the desired target state. In particular, our scheme can generate important classes of genuine multipartite entangled states, including Greenberger-Horne-Zeilinger (GHZ), Dicke and cluster states, which we discuss in detail. We show that these states can be generated with high fidelities and good heralding probabilities.

*Entanglement engine.—*We begin by describing the entanglement engine. The structure of the machine is determined by the choice of subspace, energy spectrum, and bath temperature for each qutrit, as well as the form of the interaction, all of which generally depend on the $N$-qubit target state $|\psi\rangle$. This state is obtained in a heralded manner from the steady state of the machine by projection of each qutrit to a qubit subspace. Fig. 1 shows an example targeting a GHZ state.

The machine evolution consists of a Hamiltonian contribution and a dissipative contribution due to the heat baths. The evolution is autonomous in the sense that both the Hamiltonians and the bath couplings are time independent, and the machine thus requires no work input to run. Denoting the energy basis states of qutrit $k$ by $\{|0\rangle_{k}$, $|1\rangle_{k}$, $|2\rangle_{k}\}$ and taking the corresponding energies to be $\{0,\Delta_{k}^{(1)},\Delta_{k}^{(2)}\}$, the free Hamiltonian of each qutrit is $H_{k}=\Delta_{k}^{(1)}|1\rangle_{k}\langle 1|+\Delta_{k}^{(2)}|2\rangle_{k}\langle 2|$. The free Hamiltonian of the machine is

[TABLE]

In addition, the qutrits interact via a time-independent Hamiltonian $H_{\text{int}}$, specified below.

We model the machine evolution including the heat-bath induced dissipation with a master equation of the form

[TABLE]

For simplicity, we adopt a local reset model in which the dissipator $\mathcal{L}$ corresponds to spontaneous, probabilistic, independent resets of each qutrit to a thermal state at the corresponding temperature Hartmann et al. (2006); Linden et al. (2010). That is,

[TABLE]

where $\gamma_{k}$ is the reset rate for qutrit $k$, $\tau_{k}=\exp(-H_{k}/T_{k})/\operatorname{Tr}[\exp(-H_{k}/T_{k})]$ is a thermal state of qutrit $k$, and $\otimes_{k}$ denotes tensoring at position $k$. For such a Markovian master equation description to be valid, the system-bath couplings $\gamma_{k}$ must be small relative to the system energy scale $\Delta_{k}^{(l)}$. In addition, each dissipator acts only on the corresponding qutrit, i.e. they are local. This requires that the strength of the interaction between the qutrits is at most comparable to the bath couplings $\gamma_{k}$ Hofer et al. (2017); González et al. (2017). We note that the reset model, while simple, can be mapped to a standard Lindblad-type model which can be derived from a microscopic, physical model of the baths Tavakoli et al. (2018a); Haack et al. (2019).

The goal of the machine is to produce the $N$-qubit target state by local filtering of the $N$-qutrit steady state of (2). The steady state $\rho_{\infty}$ is obtained by solving $d\rho/dt=0$, and the filter is defined by a local projection $\Pi_{k}=\openone-|R_{k}\rangle\langle R_{k}|$ of each qutrit onto the chosen qubit subspace. The state of the machine after filtering and the probability for the filtering to succeed are given by

[TABLE]

where $\Pi=\bigotimes_{k=1}^{N}\Pi_{k}$. The temperatures, filters, bath couplings $\gamma_{k}$, and the interaction must be chosen appropriately for the heralded state $\rho^{\prime}$ to approach the target state.

Here, for a given $N$-qubit target $|\psi\rangle$, we focus on the following choice for the interaction

[TABLE]

where $g>0$ is the interaction strength, and the states $|\bar{\psi}\rangle$ and $|R\rangle$ are defined by the choices of filtered qubit subspace for each qutrit. For qutrit $k$, we let $R_{k}=0,1,2$ label the level which is not part of the qubit, i.e. qubit $k$ is spanned by the two levels complementary to $|R_{k}\rangle$. Then $|\bar{\psi}\rangle$ is the embedding of the target $|\psi\rangle$ into these qubit subspaces, and $|R\rangle=|R_{1}\ldots R_{N}\rangle$. That is, $H_{\text{int}}$ swaps the target state and the state in which every qutrit is outside the filtered subspace. For example, for $N=2$, if the target state is the maximally entangled two-qubit state $|\psi\rangle=\left(|01\rangle+|10\rangle\right)/\sqrt{2}$, and we choose $|R\rangle=|20\rangle$, then the embedding into the qutrits reads $|\bar{\psi}\rangle=\left(|02\rangle+|11\rangle\right)/\sqrt{2}$.

We furhermore focus on the regime of weak inter-system coupling, where $g$ is small relative to the free energies $\Delta_{k}^{(l)}$ (where the local master equation is valid). For there to be any non-trivial evolution in this regime, the interaction needs to be energy conserving, i.e. $[H_{\text{int}},H_{\text{free}}]=0$. This restricts which target states can be generated. However, that is the only restriction. Our main result is that

Any state $|\psi\rangle$, for which the Hamiltonians $H_{\text{free}}$ and $H_{\text{int}}$ of Eqs. (1) and (5) can be constructed to satisfy $[H_{\text{int}},H_{\text{free}}]=0$, can be generated by an entanglement engine as described above.

Specifically, one may choose a single qutrit to be connected with coupling strength $\gamma_{h}$ to a hot bath at temperature $T_{h}$ and all other qubits to be connected with coupling strength $\gamma_{c}$ to cold baths at $T_{c}$. For the hot qutrit, one chooses $R_{k}=2$, while for all the cold qutrits $R_{k}=0$. The target $|\psi\rangle$ is then obtained in the limit of extremal temperatures $T_{c}=0$, $T_{h}\rightarrow\infty$, and small coupling-strength ratios $g\lesssim\gamma_{h}\ll\gamma_{c}$. A full proof is given in App. A. However, one can intuitively understand why the machine works well in this regime. When $T_{c}=0$, resets of the cold qutrits will take them to the ground state $|0\rangle_{k}$. Since for the cold qutrits $R_{k}=0$, the ground state is not part of the filtered subspace. Therefore, cold resets will only lower the filtering success probability but will not affect the overlap of the filtered state with the target state $|\psi\rangle$. Once a cold qutrit is in the ground state, the only process which can bring it back into the filtered subspace is $H_{\text{int}}$, and this can only happen once all qutrits are in the state $|R_{k}\rangle$. The hot qutrit must then be in state $|2\rangle$, which can happen via a hot reset. Hot resets also degrade the quality of the filtered state (as they destroy coherence within the filtered subspace of the hot qutrit), and hence must be much less frequent than cold reset. This way, the system is most likely to be found outside the filtered subspace (making $p_{\text{suc}}$ small), but if found inside, it is likely to be in state $|\psi\rangle$ (because it is unlikely a hot reset happens before a cold one drives the system back out). The physical intuition for the bipartite case $N=2$ was also discussed in Ref. Tavakoli et al. (2018a).

We note that, even if a given target $|\psi\rangle$ does not admit any choice of $H_{\text{free}}$ and $H_{\text{int}}$ satisfying $[H_{\text{int}},H_{\text{free}}]=0$, it may happen that by applying local unitaries to each qubit, one can obtain another state $|\psi^{\prime}\rangle$ which does. Since entanglement is preserved under local unitaries, one may then first generate $|\psi^{\prime}\rangle$ and simply apply the inverse local unitaries to obtain $|\psi\rangle$. Thus, effectively, the set of states which can be generated using the entanglement engine above consists of all states within the local unitary orbit of those $|\psi\rangle$ for which energy conservation can be satisfied.

*Energy conservation.—*We now derive conditions for $|\psi\rangle$ to admit choices of $H_{\text{free}}$ and $H_{\text{int}}$ such that $\left[H_{\text{int}},H_{\text{free}}\right]=0$. This holds if and only if every transition generated by $H_{\text{int}}$ is energy conserving w.r.t. $H_{\text{free}}$. From (5), these transitions depend on the target state and on the choice of $|R\rangle$ (which defines the filtered qubit subspaces). We can write the target $N$-qubit state as

[TABLE]

where $S_{\psi}=\{\mathbf{n}\in\{0,1\}^{N}\,|\,\langle\psi|\mathbf{n}\rangle\neq 0\}$ determines the set of basis states on which $|\psi\rangle$ has support, and $c_{\mathbf{n}}\in\mathbb{C}$. Denoting the embedding of $|\mathbf{n}\rangle$ into the $N$ qutrits by $|\bar{\mathbf{n}}\rangle$, both $|\bar{\mathbf{n}}\rangle$ and $|R\rangle$ are eigenstates of $H_{\text{free}}$ with respective eigenvalues $E_{\mathbf{\bar{n}}}$ and $E_{R}$. The conditions for energy conservation are then $E_{\mathbf{\bar{n}}}=E_{R}$ for every $\mathbf{n}\in S_{\psi}$. This can be expressed as

[TABLE]

where we have restricted to cases where the qubit states are either $\{|1\rangle_{k},|2\rangle_{k}\}$ or $\{|0\rangle_{k},|1\rangle_{k}\}$ for each qutrit (i.e. $R_{k}=0$ or $R_{k}=2$) 111Thermal resets on a given qutrit destroys entanglement with the other qutrits. For cold baths, thermal resets tend to drive the corresponding qutrit to the ground state. To suppress the effect of reset, it is therefore beneficial to choose $R_{k}=0$ when the bath temperature is cold. For infinitely hot baths, resets equalise the populations on the three levels, and it thus does not matter which subspace is filtered.. Given a target state $|\psi\rangle$, the question is thus, whether there exist choices of $R_{k}$, $\Delta^{(1)}_{k}$, and $\Delta^{(2)}_{k}$ which fulfill (Autonomous multipartite entanglement engines) for all $\mathbf{n}\in S_{\psi}$.

Although (Autonomous multipartite entanglement engines) depends only on $S_{\psi}$ and not on the coefficients $c_{\mathbf{n}}$ in (6), a general solution is not easy to obtain, because the number of variables increases with $N$. Nevertheless, (Autonomous multipartite entanglement engines) can be significantly simplified. In App. B, we show that whenever (Autonomous multipartite entanglement engines) has a solution, then it has a solution with $R_{k}=0$ for all but a single $k$. For a given $|\psi\rangle$ it is thus sufficient to check whether there exists choices of $k^{\prime}\in\{1,\ldots,N\}$, $\Delta^{(1)}_{k}$, and $\Delta^{(2)}_{k}$ fulfilling

[TABLE]

If there does, then it follows from the proof in App. A that the machine defined by these choices, with bath $k^{\prime}$ hot and all other baths cold, can generate states arbitrarily close to $|\psi\rangle$.

Below, we consider several families of genuine multipartite entangled states, important in quantum information processing, namely GHZ, Dicke and cluster states. We show that they admit solutions to (8) and hence can be generated. Furthermore, we consider the tradeoff between heralding success probability and the quality of the generated states, as well as the effect of finite temperatures, and show that they can be robustly generated also away from the ideal limit of the entanglement engine.

*GHZ states.—*We start with the GHZ state of $N$ qubits, which is commonly given as $\frac{1}{\sqrt{2}}\left(|0\ldots 0\rangle+|1\ldots 1\rangle\right)$. In this form, the state does not admit a solution to (8). However, we can instead consider $|\text{GHZ}\rangle=\frac{1}{\sqrt{2}}\left(|10\ldots 0\rangle+|01\ldots 1\rangle\right)$, which is equivalent up to a local unitary (bit flip) on the first party. One can check that $|\text{GHZ}\rangle$ does admit a solution to (8). One such solution is illustrated in Fig. 1. We take the first bath to be hot and the rest cold, and let the free Hamiltonians of the hot qutrit and each of the $N-1$ cold qutrits be

[TABLE]

To construct an energy-conserving interaction Hamiltonian, we follow the recipe above. Writing $\bar{0}$ for a string of $N-1$ zeros $0\ldots 0$, and similarly for $\bar{1}$ and $\bar{2}$, we have $|R\rangle=|2\bar{0}\rangle$. Embedding $|\text{GHZ}\rangle$ in the qutrit space, from (5) we get

[TABLE]

Once the steady state of the dynamics (2) is obtained, we apply the filter $\Pi_{h}=|0\rangle\langle 0|+|1\rangle\langle 1|$ to the hot system and the filter $\Pi_{c}=|1\rangle\langle 1|+|2\rangle\langle 2|$ to each of the cold systems. Successful filtering heralds the generation of $|\text{GHZ}\rangle$.

As explained above, the perfect GHZ state is obtained only under idealised conditions (when the temperature gradient is maximal and the coupling strength ratios tend to zero). We now consider the quality of the generated state in case of finite temperatures and varying filtering success probabilities (4). We begin with the latter.

As argued above, in the ideal limit, $\gamma_{h}\ll\gamma_{c}$, the system is most likely found outside the filtered subspace, causing $p_{\text{suc}}\rightarrow 0$ as $\gamma_{h}/\gamma_{c}\rightarrow 0$. However, away from this idealised limit, we find that the state $\rho^{\prime}$ after filtering (considered as an $N$-qubit state) may still have a high fidelity $F=\langle\text{GHZ}|\rho^{\prime}|\text{GHZ}\rangle$ with the GHZ state. Fig. 2 shows the trade-off between $F$ and $p_{\text{suc}}$ for $N=2,3,4$ systems. We see that fidelities above 90% are obtained for $p_{\text{suc}}$ at the 5%-level. Note that $p_{\text{suc}}$ is bounded, even when the fidelity is allowed to degrade. The maximal $p_{\text{suc}}$ decreases with increasing $N$, however the corresponding fidelity also increases. E.g. for $N=4$, the fidelity does not reach $F=1/2$ before $p_{\text{suc}}$ reaches its maximal value of $p_{\text{suc}}=1/9$. This suggests that as $N$ grows, the fidelity achievable up to the maximal $p_{\text{suc}}$ increases. In App. C, we derive the maximal value of $p_{\text{suc}}$ for any $N$. Finally, we note that we have also considered an analogous autonomous entanglement engine for $N=3$ with two hot systems and one cold system. However, as seen from Fig. 2, the performance in this case is worse.

We remark that, for the states considered here which have only two non-zero off-diagonal elements, a GHZ fidelity $F>1/2$ implies genuinely multipartite entanglement Gühne and Seevinck (2010). In addition, the $F>1/2$ also provides a certificate that this genuinely multipartite entanglement is strong enough to be semi-device-independently certified via the scheme of Ref. Tavakoli et al. (2018b). Furthermore, in App. D we have studied when the generated state can lead to Bell inequality violation (providing a fully device-independent certificate of entanglement).

Next, we consider the effect of finite temperatures, i.e. $T_{c}>0$ and $T_{h}<\infty$. We keep the interaction and bath coupling strengths fixed (thus also avoiding the idealised limit of vanishing couplings). The results are presented in Fig. 3. We note that even for temperatures far from the ideal limit, fidelities close to unity are possible.

Thus, our entanglement engine functions well not only in the ideal limit but also for finite temperatures and coupling strengths. In App. E, we further show that qualitatively similar results can be obtained when the simple reset model is replaced by a master equation on standard Lindblad form, which can be derived from explicit, physical modeling of the baths and interactions.

*Dicke states.—*As a second example, we consider $N$-qubit Dicke states. The Dicke state with $l$ excitations is given by

[TABLE]

where the sum is over all permutations $\sigma_{s}$ of the subsystems. Notably, setting $l=1$ returns the well-known W-states.

Again, one finds that all such states admit solutions to (8). Hence, every Dicke state can be generated by an autonomous entanglement engine. For instance, we choose the first qutrit hot and the rest cold, and the free Hamiltonians $H_{h}=\Delta^{(1)}_{h}|1\rangle\langle 1|+\Delta^{(2)}_{h}|2\rangle\langle 2|$ and $H_{c}=\Delta^{(1)}_{c}|1\rangle\langle 1|+\Delta^{(2)}_{c}|2\rangle\langle 2|$, where

[TABLE]

Note that similar solutions of (Autonomous multipartite entanglement engines) are possible also for more hot baths. For the case $(N,l)=(3,1)$, we have analytically solved the reset master equation in terms of $g,\gamma_{h},\gamma_{c}$ and computed the fidelity $F=\langle D_{1}^{3}|\rho^{\prime}|D_{1}^{3}\rangle$. Similarly, we have analytically evaluated $p_{\text{suc}}$ in (4). The tradeoff between $F$ and $p_{\text{suc}}$ is shown in Fig. 4. As for the GHZ case, we find that high fidelities can be reached with success probabilities at the few-percent level. We have also checked that increasing the number of hot systems (to two) does not improve performance.

*Cluster state.—*Finally, we consider a linear four-qubit cluster state

[TABLE]

A solution to (8) is obtained by the following free Hamiltonian, where $\Delta=\Delta^{(2)}_{c}-\Delta^{(1)}_{c}$

[TABLE]

In analogy with the previous, we consider the trade-off between the $F=\langle C|\rho^{\prime}|C\rangle$ of the generated state $\rho^{\prime}$ with the cluster state and filtering success probability $p_{\text{suc}}$. We have evaluated both $F$ and $p_{\text{suc}}$ analytically for a single hot bath, and optimised over the couplings $g,\gamma_{h},\gamma_{c}$ to obtain the results in Fig. 4. Again, high-fidelity cluster states can be generated with success probabilities at the few-percent level. Furthermore, in App. D, we have considered the device-independent certification of $\rho^{\prime}$ via Bell inequalities tailored for cluster states Scarani et al. (2005) at varying $p_{\text{suc}}$. We find that large Bell inequality violations can be obtained for every $p_{\text{suc}}$ up to its maximal value of $p_{\text{suc}}\approx 0.085$, demonstrating that the entanglement engine works well over a wide regime.

*Conclusion.—*We have given a general recipe for autonomous entanglement engines which enable heralded generation of multipartite entangled states between any number of parties. As demonstrated by several examples, a wide range of states can be targeted, including GHZ, Dicke, and cluster states. While pure target states are only generated perfectly for infinite temperature gradients and vanishing heralding success probabilities, we have explored finite temperatures and heralding probabilities as well and have found that high fidelities can be attained also away from the ideal regime.

Thus, probabilistic generation of high-quality multipartite entanglement is possible using only incoherent, thermal processes and energy-preserving interactions, requiring no work input. It would be interesting to understand if strong entanglement could be generated by an autonomous engine in a deterministic manner, i.e. without filtering. Finally, perspectives for experimental implementation could be explored. In that context, a natural question is whether genuine multipartite entangled states can be generated autonomously using only two-body Hamiltonians.

*Acknowledgements.—*We thank Marcus Huber for discussions. JBB was supported by the Independent Research Fund Denmark, AT and NB by the Swiss National Science Foundation (Grant 200021_169002 and NCCR QSIT), and GH by the Swiss National Foundation through the starting grant PRIMA PR00P2$\_$179748 .

Appendix A Autonomous generation of target states

We prove that any state $|\psi\rangle$ which admits a solution to the energy-conservation condition (Autonomous multipartite entanglement engines) can be generated by an autonomous entanglement engine. Following the main text, we write the target state as

[TABLE]

where $c_{\mathbf{n}}\in\mathbb{C}$, $\sum_{\mathbf{n}}|c_{\mathbf{n}}|^{2}=1$, and where $S_{\psi}$ is the set of binary strings $s=\{0,1\}^{N}$ such that $|\psi\rangle$ has support of $|s\rangle$. We show that the state $\rho^{\prime}$ returned by the machine described in the main text (after heralding) is indeed the target state. To this end, we must characterise $\rho^{\prime}$. For simplicity, we will first focus on the diagonal elements of $\rho^{\prime}$ and then on its off-diagonal elements.

A.1 Diagonal elements

We aim to show that the diagonal elements of $\rho^{\prime}$ correspond to the populations $|c_{\bar{\mathbf{n}}}|^{2}$, where $|\bar{\mathbf{n}}\rangle$ are the computational basis states on which the embedded target state $|\bar{\psi}\rangle$ has support. To enable the characterisation of the diagonal elements of $\rho^{\prime}$, we use flow diagrams as illustrated in Fig. 5. Such a diagram represents the transitions induced by the influence of hot and cold resets, along with the rate of said transitions, on a given support state $|\bar{\mathbf{n}}\rangle$. As illustrated; by a hot reset on $|\bar{\mathbf{n}}\rangle$ one can reach two other states, denoted by $|o\rangle$ and $|o^{\prime}\rangle$. Importantly, neither of these two states can be members of $S_{\bar{\psi}}$ since it is otherwise at odds with the conditions for an autonomous Hamiltonian. From the flow-diagram, we obtain the following steady-state condition when considering the flow into and out of the state $|o\rangle$:

[TABLE]

where we have adopted the simplified notation $P_{s}=\langle s|\rho|s\rangle$. However, since $|o\rangle,|o^{\prime}\rangle\notin S_{\bar{\psi}}$ (nor do they equal the state $|R\rangle$), they do not appear in the interaction Hamiltonian and are treated equally by the dissipation. Hence, it follows that $P_{o}=P_{o^{\prime}}$. This leads us to re-write (19) as

[TABLE]

Let us now consider the filtered subspace, i.e. the space in which the heralded state $\rho^{\prime}$ lives. Since the filtering corresponds to projecting each qutrit onto a qubit subspace, there are consequently $2^{N}$ computational basis states spanning the filtered sub-space. Of these, $\nu=|S_{\bar{\psi}}|$ are members of $S_{\bar{\psi}}$, whereas another $\nu$ are reachable by a hot reset to each element in $S_{\bar{\psi}}$. Denote the latter set of states by $G_{h}$. The remaining $2^{N}-2\nu$ states have no population (diagonal element equal zero) since they can neither be reached via the interaction Hamiltonian nor via resets. Let $\bar{P}_{o}$ denote renormalised $P_{o}$ after filtering, i.e., $\bar{P}_{o}=\langle o|\rho^{\prime}|o\rangle$. Normalisation requires that

[TABLE]

However, due to the symmetries of the interaction Hamiltonian and the linearity of the dynamics, we may write $\bar{P}_{o}=|c_{o}|^{2}\bar{P}_{\text{S}}$ for $o\in S_{\bar{\psi}}$ for some constant population $\bar{P}_{\text{S}}$ independent of $o$. Similarly, we may write $\bar{P}_{o}=|c_{o}|^{2}\bar{P}_{\text{G}}$ for $o\in G_{h}$ for some constant population $\bar{P}_{\text{G}}$ independent of $o$. The normalisation condition reduces to

[TABLE]

which together with (20) gives

[TABLE]

In the limit $\gamma_{h}\ll\gamma_{c}$ we have $\bar{P}_{\text{S}}\rightarrow 1$, and therefore also $\bar{P}_{\text{G}}\rightarrow 0$. Consequently, we have found that in the given limit, for $\bar{\mathbf{n}}\in S_{\bar{\psi}}$

[TABLE]

These are the desired diagonal elements.

A.2 Off-diagonal elements

We now aim to show that the off-diagonal elements of $\rho^{\prime}$ correspond to $c_{n}c_{n}^{*}$. Due to hermiticity, it is sufficient to consider the upper triangle in the matrix of $\rho^{\prime}$. Among these off-diagonal entries, there are $\binom{\nu}{2}$ that correspond to coherences generated between the computational basis states associated to $n,n^{\prime}\in S_{\bar{\psi}}$ (we have dropped the notation in bold ($\bar{\mathbf{n}}$) since in this section $n$ will sometimes be a member of $S_{\bar{\psi}}$). Another $\nu$ off-diagonals correspond to coherences generated between the computational basis states assciated to $n\in S_{\bar{\psi}}$ and the state $|R\rangle$. The remaining off-diagonal elements are not reachable by the dynamics (neither via resets nor via the Hamiltonian) and therefore equal zero. We use the short-hand notation $\rho_{n,n^{\prime}}=\langle n|\rho|n^{\prime}\rangle$ to write the reset master equation in the steady state as

[TABLE]

For the first term in Eq. (25) we have that

[TABLE]

Taking $n,n^{\prime}\neq R$, the two middle terms vanish. Moreover, if $n,n^{\prime}\notin S_{\bar{\psi}}$ also the first and fourth term vanish. If $n,n^{\prime}\in S_{\bar{\psi}}$ then we have $\langle n|\bar{\psi}\rangle=c_{n}$ and $\langle\bar{\psi}|n^{\prime}\rangle=c_{n^{\prime}}^{*}$ and therefore $\langle n|[H,\rho]|n^{\prime}\rangle=g\left(c_{n}\rho_{R,n^{\prime}}-c_{n^{\prime}}^{*}\rho_{n,R}\right)$. Thus,

[TABLE]

For the second term in Eq. (25) a direct calculation gives

[TABLE]

where the bar-sign denotes $\bar{s}=s_{2}\ldots s_{N}$. Moreover, the third term in (25) straightforwardly evaluates to

[TABLE]

where $\overleftarrow{s}=s_{1}\ldots s_{k-1}$ and $\overrightarrow{s}=s_{k+1}\ldots s_{N}$. Notice that this term vanishes for $k=2,\ldots,N$ if either $n$ or $n^{\prime}$ are members of $S_{\bar{\psi}}$. In conclusion, for $n,n^{\prime}\in S_{\bar{\psi}}$, we can re-write (25) as

[TABLE]

When $n\neq n^{\prime}$ (since one cannot transition between two support states by a hot reset) Eq. (28) becomes $\delta_{n_{1},n_{1}^{\prime}}\sum_{j}\rho_{j\bar{n},j\bar{n}^{\prime}}=\delta_{n_{1},n_{1}^{\prime}}\rho_{n,n^{\prime}}$. Furthermore, by hermiticity we have that $\rho_{R,n^{\prime}}=\rho_{n^{\prime},R}^{*}$, and due to the symmetries of the Hamiltonian it also holds that $\rho_{n,R}=c_{n}L$ where $L$ is a constant related to the population in the steady-state that is independent of $n$. With this in hand, we consider the three equations obtained from (30):

[TABLE]

where in the first equation we have taken $n,n^{\prime}\in S_{\bar{\psi}}$ with $n\neq n^{\prime}$, in the second equation we have taken $n,n^{\prime}\in S_{\bar{\psi}}$ with $n=n^{\prime}$, and in the third equation we have taken $n,n^{\prime}\in S_{\bar{\psi}}$ with $n=n^{\prime}$ but then replaced $n_{1}$ with the index $i$ which runs over the two values $i\neq n_{1}$. Summing over $i$ in the equation (33) gives

[TABLE]

Inserted into the equation (32) we obtain

[TABLE]

Finally, when inserted into the equation (31), we can obtain the off-diagonal elements from the diagonal elements of $\rho^{\prime}$. We obtain

[TABLE]

However, the ratios between the off-diagonal terms are conserved after filtering if they belong to the filtered subspace. We use the notation $\bar{\rho}_{s,s^{\prime}}=\langle s|\rho^{\prime}|s^{\prime}\rangle$. Then, taking the relevant limit of $\gamma_{h}\ll\gamma_{c}$, we obtain

[TABLE]

The right-hand-side features a diagonal element which was evaluated in (24). In the relevant limit, we obtain the final result

[TABLE]

In conclusion, we have shown that the heralded state $\rho^{\prime}$ is the target state.

Appendix B Simplified conditions for energy conservation

B.1 A single hot system is sufficient

Here, we show that if the conditions (Autonomous multipartite entanglement engines) for the interaction to be energy conserving can be solved using $q$ hot systems (i.e. systems with $R_{k}=2$) and $N-q$ cold systems (i.e. systems with $R_{k}=0$), then there also exists a solution with just a single hot system and $N-1$ cold systems.

To prove this, we show that any set of valid energies $\Delta_{k}^{(1)}$, $\Delta_{k}^{(2)}$ fulfilling the energy-conservation condition for $q$ hot systems allows one to define another set of energies $\{\varepsilon_{k}^{(1)}$, $\varepsilon_{k}^{(2)}\}$ which fulfill the corresponding condition with a single hot system. Without loss of generality (as one may always permute the parties), we can take the hot systems to be the first ones. Then the energy-conservation condition with $q$ hot systems reads

[TABLE]

while the corresponding condition with a single hot system ($q=1$) becomes

[TABLE]

Note that the energies must satisfy $\Delta_{k}^{(2)}>\Delta_{k}^{(1)}>0$ and similarly $\varepsilon_{k}^{(2)}>\varepsilon_{k}^{(1)}>0$. To construct a solution to (40) given a solution to (39), we choose

[TABLE]

and

[TABLE]

for some $t_{k}$ satisfying $t_{k}>\Delta_{k}^{(2)}$. Note that with these choices we have $\varepsilon_{k}^{(2)}>\varepsilon_{k}^{(1)}>0$ for $k=2,\ldots,N$, as desired. Inserting in (40), we get

[TABLE]

This reduces to (39) provided that

[TABLE]

which is solved by

[TABLE]

It is easy to see that $\varepsilon_{1}^{(2)}>\varepsilon_{1}^{(1)}>0$. We thus have a valid choice of energies $\varepsilon_{k}^{(1)}$, $\varepsilon_{k}^{(1)}$ for which (40) reduces (39). Hence, any solution with $q$ hot systems also implies the existence of a solution with a single hot system, as claimed.

B.2 Identical energy structures for all hot and all cold systems

If the energy spectra of all hot systems (i.e. all systems with $R_{k}=2$) are identical, and similarly those of cold systems (with $R_{k}=0$) are identical, then the energy-conservation conditions can be simplified. Note that all the examples given in the main text (for GHZ, Dicke, and cluster states) belong to this setting.

Specifically, here we show that if $\Delta^{(1)}_{k}$ and $\Delta^{(2)}_{k}$ depend only on $R_{k}$, then the existence of $R\in\{0,2\}^{N}$ and a choice of energies fulfilling (Autonomous multipartite entanglement engines) in the main text is equivalent to the existence of a vector $\mathbf{r}\in\{0,1\}^{N}$ such that $\mathbf{r}\neq\mathbf{0},\mathbf{1}$ and for each pair of vectors $\mathbf{n},\mathbf{n}^{\prime}\in S_{\psi}$ either

[TABLE]

or

[TABLE]

where $c<0$ is a constant independent of $\mathbf{n},\mathbf{n}^{\prime}$, and $\mathbf{0}=(0,0,\ldots,0)$ and $\mathbf{1}=(1,1,\ldots,1)$. That is, the interaction (5) can be made energy conserving if and only if an $\mathbf{r}\neq\mathbf{0},\mathbf{1}$ exists fulfilling (53)-(54).

Before we proceed with the proof, we illustrate (53)-(54) and the notation introduced above in the simplest setting of two parties. We take the maximally entangled state $|\Psi^{+}\rangle=\left(|01\rangle+|10\rangle\right)/\sqrt{2}$ as the target and choose $|R\rangle=|20\rangle$. The target has support on just two states, $S_{\psi}=\{\mathbf{n},\mathbf{n}^{\prime}\}$, where

[TABLE]

It is straightforward to verify that (54) is satisfied for $\mathbf{r}=(1,0)$, with $c=-1$. Hence, $|\Psi^{+}\rangle$ can indeed be generated autonomously. Looking at (Autonomous multipartite entanglement engines), we see that the conditions on the energies coming from $\mathbf{n}$ and $\mathbf{n}^{\prime}$ are respectively

[TABLE]

and

[TABLE]

Thus, the two qutrits have the same maximal energy but inverted level structures. The gap between the two lower levels for the second qutrit equals the gap between the upper two levels for the first qutrit. This corresponds exactly to the entanglement engine of Ref. Tavakoli et al. (2018a).

The conditions (53)-(54) can be defined as follows. If we define a vector $\mathbf{r}\in\{0,1\}^{N}$ such that $r_{k}=0$ if $R_{k}=0$ and $r_{k}=1$ for $R_{k}=2$, then for each $\mathbf{n}\in S_{\psi}$, the condition $E_{\mathbf{n}}=E_{\bar{R}}$ from (Autonomous multipartite entanglement engines) of the main text can be expressed as

[TABLE]

The question is, whether there exist choices of $\mathbf{r}$, $\Delta^{(1)}_{k}$, and $\Delta^{(2)}_{k}$ which fulfill this. Rewriting, we have

[TABLE]

Now, if the energy structures of all qutrits with the same $R_{k}$ are the same, then the energies appearing under each sum become independent of $k$. Let us denote the energy gaps of qutrits with $R_{k}=0$ by $\delta_{1}=\Delta^{(1)}_{k}$ and $\delta_{2}=\Delta^{(2)}_{k}-\Delta^{(1)}_{k}$ and those of qutrits with $R_{k}=2$ by $\delta_{3}=\Delta^{(1)}_{k}$ and $\delta_{4}=\Delta^{(2)}_{k}-\Delta^{(1)}_{k}$. Then (59) becomes

[TABLE]

which is equivalent to

[TABLE]

where $\mathbf{1}=(1,\ldots,1)$ and $|\mathbf{r}|$ is the number of 1’s in $\mathbf{r}$. This must hold for every $\mathbf{n}\in S_{\psi}$, and thus we have a set of linear equations

[TABLE]

where $\nu$ is the number of elements of $S_{\psi}$. Regarding $\mathbf{\delta}=(\delta_{1},\ldots,\delta_{4})$ as a variable, we would like to know when there exists $\mathbf{r}\in\{0,1\}^{N}$ such that (62) has a solution over $(\mathbb{R}^{+})^{4}$, i.e. a positive solution. Given such a solution, for any $l=1,\ldots,\nu$ we must have

[TABLE]

where the right-hand side is independent of $l$. Note that this condition can never be fulfilled if $\mathbf{r}=\mathbf{0}$ or $\mathbf{r}=\mathbf{1}$, because the two sides of the equation then have opposite signs. However, if the condition is satisfied, then for any pair $l,l^{\prime}=1,\ldots,\nu$ we have

[TABLE]

Hence, for a positive solution to exist, for each pair of support states either $\mathbf{n}^{(l)}$ and $\mathbf{n}^{(l^{\prime})}$ have an equal number of 1’s in positions where $\mathbf{r}$ has 0 and an equal number of 1’s in positions where $\mathbf{r}$ has 1, or

[TABLE]

is a negative constant independent of $l$, $l^{\prime}$. On the other hand, if an $\mathbf{r}\neq\mathbf{0},\mathbf{1}$ exists fulfilling these conditions, then a positive solution of (62) is guaranteed to exist. This is because the left-hand side of (63) is then independent of $l$ and thus one can always find positive $\delta_{1}$ and $\delta_{4}$ which make the equality true.

Appendix C Maximal filtering probability in GHZ-state machine

Naturally, since $N$ local filters are performed on the steady state of an $N$-qutrit autonomous thermal machine, the probability of a successful filtering decreases with $N$. It is therefore reasonable to ask what this maximal possible success probability is. This can be determined analytically by considering the flow of population in the steady state of the GHZ machine.

Since a cold reset always takes a system out of the filtered sub-space, the maximal success probability is obtained in the limit $\gamma_{h}\gg\gamma_{c}$, i.e. the opposite of the limit maximising the fidelity of the generated state with the target state. To determine $p_{\text{suc}}$ in this limit, let $S_{k}$ denote the set of all eigenstates of the joint free Hamiltonian where $k$ cold qutrits are in one of the excited states (all in the same one), while the remaining $N-k-1$ cold qutrits are in the ground state. For instance, in $S_{N-1}$ we have the states $S_{N-1}=\{|0,\bar{1}\rangle,|1,\bar{1}\rangle,|2,\bar{1}\rangle,|0,\bar{2}\rangle,|1,\bar{2}\rangle,|2,\bar{2}\rangle\}$ while $S_{0}$ consists of the states $S_{0}=\{|0,\bar{0}\rangle,|1,\bar{0}\rangle,|2,\bar{0}\rangle\}$. We will compare the flows of population into and out of the $S_{k}$. However, first we argue that within each $S_{k}$, the populations on each of the states are equal in the steady state. We first note that all processes (the evolution driven by the $H_{\text{int}}$ of the GHZ machine, as well as hot and cold resets) are symmetric in the states $|\bar{1}\rangle$ and $|\bar{2}\rangle$ of the cold qutrits. The populations of states with the hot qutrit in a fixed state and a fixed number of cold qutrits excited to the same excited state, and which differ only in whether this state is $|1\rangle$ or $|2\rangle$ must therefore be equal in the steady state. In contrast, $H_{\text{int}}$ is not symmetric in the states $|0\rangle$, $|1\rangle$, $|2\rangle$ of the hot qutrit, and hence populations of states with the hot qutrit in different levels are not expected to be equal in the steady state in general. However, in the limit $\gamma_{h}\gg\gamma_{c}$, there are many hot resets between each cold one. This will then equalise the populations within each set $S_{k}$ before a cold reset causes a transition to $S_{k-1}$. Hence, all populations with each $S_{k}$ are equal in the steady state.

We can now draw the flow diagram shown in Fig. 6 for population transfer between the $S_{k}$. In the steady state, the flow into each set $S_{k}$ must equal the flow out. If we denote the population per state in $S_{k}$ by $P_{k}$, we therefore have, for $k=2,\ldots,N-1$

[TABLE]

The number of states in the set $S_{k}$ is given by

[TABLE]

Inserting in (66) and rearranging, one finds that

[TABLE]

From this it follows that $P_{N-1}=P_{1}$ and $P_{N-2}=P_{2}$ etc. That is,

[TABLE]

To determine the relation with $P_{0}$, we note that $H_{\text{int}}$ drives swaps between the states $|2\bar{0}\rangle\leftrightarrow\frac{1}{\sqrt{2}}(|1\bar{1}\rangle+|0\bar{2}\rangle)$ and hence between $S_{0}$ and $S_{N-1}$. This process is a unitary rotation. Nevertheless, in the steady state it still results in a flow of population with a constant rate, which we can denote $\nu_{g}$. Focusing on the flow in and out of $S_{N-1}$, we can write

[TABLE]

As argued above, when $\gamma_{h}\gg\gamma_{c}$, all states in each $S_{k}$ are equally probable, and so

[TABLE]

Now, if further $\nu_{g}\gg(N-1)\gamma_{c}$, then

[TABLE]

Finally, normalisation of the steady state requires that

[TABLE]

Together, Eqns. (68), (72), and (73) provide $N$ independent equations from which the populations $P_{k}$, $k=0,\ldots,N-1$ can be determined. Explicitly, we can first express everything in terms of $P_{0}$. For $k\geq 1$

[TABLE]

where we used that $P_{1}=P_{N-1}=P_{0}$. Then, from (73)

[TABLE]

and hence

[TABLE]

where $h_{n}$ is the $n$’th harmonic number. We can now compute the probability for successful filtering, given the steady-state populations (74) and (78). The success probability becomes

[TABLE]

where the last line is valid for large $N$. We note that the assumption $\nu_{g}\gg(N-1)\gamma_{c}$ leading to (72) may not formally be justified for the local master equation. However, we have checked that the final expression (80) is consistent with solutions obtained for $N\leq 8$ without making this assumption.

It is interesting to observe that the critical $p_{\text{suc}}$ for obtaining a non-trivial GHZ-state fidelity approaches the above maximal value (79) of $p_{\text{suc}}$ rapidly already for $N=3,4$ displayed in Fig. 2. Provided that this observation extends to larger $N$, it is interesting to note that genuinely multipartite entanglement can be generated with a success probability which decreases only log-linearly with $N$.

Appendix D Nonlocality versus filtering probability in the GHZ-state and cluster-state machines

A particularly strong form of entanglement is that which can violate a Bell inequality. Therefore, we have considered whether the states generated by the GHZ machine at fixed success probabilities have the ability of violating Bell inequalities. To this end, we have focused on the Mermin inequalities Mermin (1990) which is a family of Bell inequalities applicable to scenarios in which $N$ observers share a state and perform one of two local measurements with binary outcomes. These inequalities are known to be maximally violated by a GHZ state. Let the input of the $k$’th observer in the Bell scenario be $x_{k}\in\{0,1\}$ and the corresponding output be $a_{k}\in\{0,1\}$. We use a somewhat modified variant Tavakoli (2016) of the Mermin inequalities which reads

[TABLE]

where

[TABLE]

We have fixed the measurements of each observer to be those required for a maximal violation with a GHZ state. For $N=2$, the optimal measurements are $\sigma_{x}$ and $\sigma_{z}$ for one observer, and $(\sigma_{z}+\sigma_{x})/\sqrt{2}$ and $(\sigma_{z}-\sigma_{x})/\sqrt{2}$ for the other observer. For $N=3$ we have let all three observers perform either $\sigma_{x}$ or $\sigma_{y}$, and for $N=4$ one observer performs either $\sigma_{x}$ or $\sigma_{y}$ whereas the remaining three choose between $(\sigma_{x}+\sigma_{y})/\sqrt{2}$ and $(\sigma_{x}-\sigma_{y})/\sqrt{2}$. We have numerically obtained the trade-off between nonlocality the filtering success probability. The results are illustrated in Fig. 8. We conclude that the states generated by the GHZ machine can violate Bell inequalities for reasonable $p_{\text{suc}}$.

We have also performed an analogous analysis for the states generated at fixed success probabilities in the cluster-state machine. Specifically, we have considered whether these states can violate a Bell inequality tailored for cluster states Scarani et al. (2005). We have restricted ourselves to the measurements optimal for a cluster state $1/2(|0000\rangle+|0011\rangle+|1100\rangle-|1111\rangle)$ which is unitarily equivalent to the target state (15). Hence, after a suitable local unitary, the Bell expression reads

[TABLE]

which is bounded by $B\leq 2$ in all local hidden variable models. With a cluster state, one can achieve $B=4$. The trade-off between $B$ and the success probability of fitering is displayed in Fig. 8. We find that the generated states are nonlocal for any $p_{\text{suc}}$ up to its maximal value.

Appendix E Lindblad-type master equation

To demonstrate that our results are not restricted to the simple reset model employed in the main text, here we provide a Lindblad-type master equation, which can be derived from a microscopic model with bosonic baths. The reset model (2)-(3) is replaced by

[TABLE]

where $\Gamma_{k}$ denotes the rate of a transition, $n_{B}(E,T)=1/(e^{E/T}-1)$ is the Bose-Einstein distribution and $\mathcal{D}$ denotes the dissipator 222We remark that this model suppresses one transition, in analogy with the Lindblad-type master equation considered in Ref. Tavakoli et al. (2018a).

Results from the Lindblad-type model qualitatively agree with those of the reset model. As an example, we again consider a GHZ target state for three parties ($N=3$), solve for the steady state, and find the GHZ fidelity of the filtered state as a function of $T_{h}$ and $T_{c}$. The result is shown in Fig. 9. Just as in Fig. 3 in the main text, we see that high fidelities can be attained with reasonably low temperature gradients. Parameter values are chosen based on recent experimental results in circuit-QED Pop et al. (2014); Jerger et al. (2016); Cottet et al. (2017); Lin et al. (2018), see also Tavakoli et al. (2018a).