Entanglement and its dynamics in open, dissipative systems

Quantum mechanical entanglement can exist in noisy open quantum systems at high temperature. A simple mechanism, where system particles are randomly reset to some standard initial state, can counteract the deteriorating effect of decoherence, resulting in an entangled steady state far from thermodynamical equilibrium. We present models for both gas-type systems and for strongly coupled systems. We point out in which way the entanglement resulting from such a reset mechanism is different from the entanglement that one can find in thermal states. We develop master equations to describe the system and its interaction with an environment, study toy models with two particles (qubits), where the master equation can often be solved analytically, and finally examine larger systems with possibly fluctuating particle numbers. We find that in gas-type systems, the reset mechanism can produce an entangled steady state for an arbitrary temperature of the environment, while this is not true in strongly coupled systems. But even then, the temperature range where one can find entangled steady states is typically much higher with the reset mechanism.


Introduction
In quantum information theory entanglement between parts of a system has been identified as the key resource that can possibly make quantum information processing more powerful than classical information processing. Entanglement can also be a resource for long-distance quantum communication or distributed quantum computation, and it is at the heart of some quantum communication protocols. But entanglement is fragile under the influence of environment induced decoherence. All engineering hence thrives to better control and manipulate the quantum information stored in the system while keeping the detrimental effects of decoherence low.
In nature, on the other hand, we mostly find less controllable systems, especially if the system size becomes macroscopic as in gases, fluids, solids or even biological systems. Since these systems are usually open, noisy systems at possibly high temperatures one expects that environment induced decoherence will erase all entanglement between system degrees of freedom. This reasoning is true except for three cases.
First, the environment and its coupling to the system could be special in a way that it creates rather than destroys entanglement. However, it is unlikely to find such an environment in nature where usually thermalization dominates, and we will only briefly touch upon the subject of such environments in this paper.
Second, if the system has an entangled ground state, as many solid state systems do, its thermal state will be entangled in a certain temperature range above zero by a continuity argument. Coupling to a heat bath drives a system into its thermal state. But there is a temperature threshold for the bath above which the thermal state will be unentangled. Third, the system might have a built-in entropy drain, meaning that the correlations with the environment are, by one way or another, erased such that the system can re-build entanglement through its quantum mechanical interactions. This entropy drain may even be local to exclude the trivial cases where entanglement is simply "pumped" into the system, e.g., by injecting fresh, entangled Bell pairs.
In [1] we proposed such a local entropy drain in form of a reset mechanism, where system particles are randomly replaced by particles in some standard, mixed state of sufficiently low entropy. Note that such a mechanism cannot create entanglement, on the contrary, it erases any entanglement that might still be present between the particle that is reset and the rest of the system. Only the interplay with the system Hamiltonian can lead to entanglement in a steady state that is possibly far from thermodynamic equilibrium. This reset mechanism was studied for a toy model with two qubits, where analytic solutions could be obtained. Also a multipartite scenario for a (simplified) gas model was discussed, and further generalizations were suggested. By gas-type systems we mean systems in which the decoherence processes act locally on the system particles, by strongly coupled systems we mean those where the decoherence processes act globally. To be more precise, local decoherence processes are those, which induce transitions between the eigenstates of the local, free Hamiltonian alone, while global decoherence processes induce transitions between eigenstates of the total Hamiltonian.
In this paper we review the key idea of a reset mechanism but provide more in-depth material than in [1]. We elaborate on the generalizations suggested in [1], namely on the influence of local entropy drains on the dynamics of entanglement and on the steady-state entanglement in gas-type systems as well as in strongly coupled systems.
We prove that the master equation describing the evolution of the system coupled to a heat bath and subject to a reset mechanism is of Lindblad form and hence generates a completely positive, i.e., physical map. We analytically solve the master equation for small systems of two spins with special interaction Hamiltonians, which enables us to illustrate the main features of the reset mechanism. In particular we show the following.
(i) Steady-state entanglement in systems with reset mechanism is different from the entanglement in thermal states. (ii) In strongly coupled systems with constant coupling steady-state entanglement with reset can exist for higher temperatures than the entanglement in the thermal state, which is the steady state without reset. (iii) In gas-type systems steady-state entanglement with reset can exist even for arbitrary temperatures.
These features are not due to the specially chosen interaction Hamiltonians and decoherence processes. We demonstrate that the above properties are almost independent of both. One can also relax the conditions on the reset states and take mixed states with sufficiently low entropy instead of pure states. Finally, a generalization to larger system sizes, possibly even with fluctuating particle numbers, still leads to similar results. Hence, the reset mechanism is at the same time simple and generic.
We remark that in cavity QED an incoherent generation of entanglement has been proposed, which bears resemblance to our work [2]. There, an atom couples to two leaky optical cavities and is driven by a white noise field. This incoherent driving can, when the atom is finally traced out, result in entanglement between the cavity modes. Entanglement is generated for intermediate cavity damping rates and intensities of the noise field, an effect labeled "stochastic resonance" in [2,3]. We believe that this effect is more correctly interpreted as an example for a reset mechanism. In a subsequent work [4] strongly related to [2], one single cavity entangles two atoms, giving yet another example for a reset mechanism even closer to the setups of this paper.
The reset mechanism is certainly not a preferred way to actively protect entanglement and mostly cannot even be compared to such strategies, but, because of its simplicity and generality, there is hope that such a mechanism may ultimately be identified in natural processes leading to an increased understanding whether entanglement can play a role in systems at high temperatures.
The paper is organized as follows. We first concentrate on simple models with only two particles (qubits). In section 2 we motivate the description by a master equation, explain in which cases the model is valid, and study several specific Hamiltonians and noise channels analytically and others numerically. We also compare entangled steady states resulting from a reset mechanism to entangled steady states resulting from special choices of interaction Hamiltonian and decoherence process. We show in section 3 that we can find the same features in strongly coupled systems, and we give the conditions to be met by the reset mechanism such that entangled steady states can exist. Then, in section 4, we extend the model to include more qubits and discuss the meaning of different kinds of entanglement that we use. Finally we give a summary of the results in section 5.

Gas-type systems
In this chapter we discuss a toy model with only two particles, which we take as spin-1/2 systems or qubits for simplicity. The toy model shows all the features that we will later find in larger systems and it has the advantage that we can show many results analytically leading to an increased understanding of the involved processes. We will formulate the equations for an arbitrary number N of qubits, so that we can refer to them later in section 4.

Master equation for gas-type systems
In a gas particles are weakly coupled in the sense that most of the time they do not considerably interact with each other unless they collide. In the meantime they only feel their local, free Hamiltonian and are subject to individual, local decoherence processes, e.g. through interactions with thermal photons (radiative damping). If we pick a subset consisting of N gas particles and consider these as the system, collisions with the remaining gas particles are another source of decoherence (nonradiative damping or dephasing). In a master equation that models this gas-type system we replace the original, time-dependent collision Hamiltonian by an averaged, time-independent interaction Hamiltonian. Since the interaction Hamiltonian does not modify the energy landscape in this model, the local, radiative decoherence processes tend to drive the system to the thermal state of the free Hamiltonian, for which we choose the form We leave the interaction Hamiltonian H unspecified for the moment. For two qubits, we will often use the Ising Hamiltonian for analytic discussions, whereas more complicated Hamiltonians will be treated numerically. We write the total Hamiltonian as H total = H + H free such that the master equation iṡ where L noise is a Liouville operator representing the noise channels. We describe the noise channels by the Lindblad operator [5,6] where σ ± = (σ x ± iσ y )/2 and the σs are Pauli operators. Parameters B and C give the decay rate of inversion 1+σz 2 and polarization σ ± under the action of L noise , and s = lim t→∞ (1 + σ z )/2 t = (e ωβ + 1) −1 ∈ [0, 1] depends on temperature, where s = 1/2 corresponds to T = 1/β = ∞ (we set the Boltzmann and Planck constant equal to one). The definition of s stems from laser physics where inversion occurs corresponding to "negative temperatures". Many authors usen = 1/(e ωβ − 1) instead of s andn+1 instead of (1−s). Then, no negative temperatures are possible. The noise channel is derived assuming certain approximations, e.g. the Markov approximation. Note, however, that this is not an essential assumption as we will demonstrate later in an example (see Figure 2 and related text).
An important special case of (4), obtained by setting B = 0 and C = 2γ, is the dephasing channel well-known especially in its integrated form as a completely positive map As with the Ising Hamiltonian (2), we will often use the dephasing channel for analytic discussions because of its simplicity.
We are interested in the steady state of this master equation for N = 2 qubits and the question whether there is entanglement in this state. At this point we simply state the following results since we will later solve a more general master equation that contains equation (3) as a special case.
As an easy example we start with the Ising interaction Hamiltonian (2). The steady state will be the tensor product of the thermal states of each free Hamiltonian H z since this state commutes with H Ising and since dephasing noise does not change the diagonal elements of the density matrix. In conclusion, we have the unentangled steady state When can we hope to find an entangled steady state? As we know the radiative decoherence processes drive the system into the thermal states of the free Hamiltonian. If the interaction Hamiltonian can entangle these states we may find an entangled steady state at least for low temperatures. As an example, consider an interaction Hamiltonian H = gσ (1) x σ (2) x while the other terms stay the same as in the example above. We set C = 1/2 B (no non-dissipative processes). The steady-state density matrix is then This density matrix can be entangled. We measure the entanglement between two sets of qubits by the negativity [7,8], which is given with respect to a bipartition A-Ā as N A = (||ρ TA || 1 − 1)/2, where T A means the partial transpose with respect to A. For two qubits, we omit the label A since there is only one bipartition, and the negativity can assume values between zero (separable state) and 1/2 (maximally entangled state). The reason why we choose the negativity as a measure throughout the paper is that we will use a generalization thereof in the multipartite case where the generalization of other entanglement measures might be hard to compute.
The negativity of the state above is which can be larger than zero but will always vanish for high temperatures of the bath, s → 1/2 (see Figure 1). We see that the steady state of the master equation (3) can be entangled for specially chosen interaction Hamiltonians, but only below a certain temperature threshold. Furthermore, if the free, local Hamiltonian is too strong, i.e., ω dominates by far all other parameters, there is also no entanglement. This statement applies to other models involving the Hamiltonian (1) as well, and, accordingly, ω should have the same order of magnitude as the other parameters.

Example for a gas-type system: the spin gas
Spin gases [9,10] are an example for such gas-type systems. A spin gas is a system of quantum spins with stochastic, time-dependent interactions. A physical model of a spin gas is a system of N classically moving particles with additional, internal and free Hamiltonian The decoherence processes are given by equation (4). We choose the decay rate of inversion, B, as inverse unit timescale and the other parameters as C = B/2 (no dephasing noise) and ω = B. The parameter g on one of the axes is also measured in units of B, whereas s is dimensionless.
spin degrees of freedom. Upon collision, these quantum degrees of freedom interact according to some specified Hamiltonian. In [9,10] the interaction Hamiltonians were chosen locally unitarily equivalent to the Ising interaction leading to a description in terms of weighted graph states. Hence, in such spin gases, classical kinematics drives the evolution of the quantum state, and also the decoherence of arbitrary probe systems put into the gas and subjected to interactions with it. In general, multiple non-consecutive collisions of particles are possible. The spin gas remembers its whole interaction history, and it provides a microscopic model with non-Markovian decoherence.
Assume that we have two selected gas particles (e.g. another species) that we consider as the system, while the other gas particles act as the environment that induces decoherence when interacting with the system particles. The rare interactions between the two system particles happen only during the short times when they collide. In the longer times in between they are not coupled and subject to local decoherence processes, i.e., interactions with the environment gas particles. The induced decoherence processes are equivalent to dephasing channels (corresponding to B = 0, C = 2γ in equation (4)). In such a situation, any entanglement between the two qubits that may either have been present initially or have built up on a short time scale will eventually be destroyed by the interactions with the other gas particles [1].

The spin gas with reset mechanism
For the moment, we stick to the toy-model with only to selected gas particles. Imagine now that the two particles can, at a certain rate, leave the box in which the gas is contained and are instantly replaced by fresh qubits that are in a standard mixed state with sufficiently low entropy. Instead of a replacement of system particles one can equivalently picture a measurement of the particle and a subsequent preparation in this standard state. Note that the last step need not be an active procedure but can, e.g., result from a spontaneous decay to this state. We call both procedures a reset mechanism. Certainly, by a reset, we did not introduce entanglement into the system always consisting of two qubits. On the contrary, any entanglement that might have been present between the particle that has left the box and the one that is still inside leads to a description of the latter by a more mixed density matrix (closer to the identity). But the advantage is that we have lowered the local entropy of the system since the new particle has no correlations with the environment. This new particle can then become entangled with the other one on short time scales. We said above that in a spin gas with zero rate of qubit exchange the steady state will not be entangled. For infinite exchange rate the system would always be in a pure (or a standard separable) state and there is also no entanglement. If, however, the rate at which the qubits leave the spin gas is in a certain intermediate parameter regime one can hope that there is entanglement in the system on average. Here, averaging means taking the mean density matrix of many simulation runs. Later, the solutions to master equations are assumed to resemble the evolution of such a mean density matrix and, for explanations of certain (entanglement) features, this picture will sometimes be invoked. Note that there are also other ansatzes. In [11] the solution of the master equation represents a smoothed version of a single simulation run, where smoothing is achieved by a timeintegration kernel. The solution of the master equation does then not follow the rapid changes of the single density matrix, but sees only the slower changes resolved by a so-called coarse-grained timescale, which is related to the support of the integration kernel. Figure 2 indeed shows entanglement in steady states in a simulation of a spin lattice gas with an Ising-type interaction [9,10].
The above scenario with two qubits might seem a little artificial. However, if we extend the setup to systems with more qubits and allow fluctuating particle numbers we can drop the requirement that selected particles must be instantaneously replaced. We can regard the spin-gas particles that do not belong to the system as a "hot bath" introducing decoherence and destroying entanglement in the system, and the source of qubits in standard states as a "cold bath" that can counteract the effect of the hot bath and preserve entanglement in the steady state. The analogy with a cold bath has some limitations as we will point out in section 3. We will deal with multipartite qubit systems in section 4. For the moment, we will stick to the two qubit system, for which we can find a master equation that reflects the properties of the "cold bath" and that we can solve analytically in certain cases. Observe that the master equation will again incorporate the Markov-assumption, whereas the spin gases are non-Markovian systems [9,10] and also partly have non-local decoherence processes. The essential features on the other hand will be qualitatively the same in some parameter regimes of the spin gas, where these effects play a minor role.
In the next subsection, our goal is to transfer the idea of a reset mechanism from the specific example of a spin gas to a description in form of a master equation, suitable for any spin system. . Steady-state entanglement between two selected qubits as measured by the negativity in a spin lattice gas (8×8 lattice, 2+18 qubits). The probability that a particle is exchanged for a fresh one in one time step of the simulation is plotted on the horizontal axes. Hence, the value 1 corresponds to an infinite exchange rate. The special qubits interact 1000 times stronger with each other than with the 18 qubits that form the environment, i.e., as physical particles, they are e.g. of a different kind than the environment particles. The density matrices from which the negativity is derived were averaged over 10000 simulation runs. Other parameters: initial distance between special qubits: 1, interaction phase picked up during a collision between them: ψ = 0.1; interaction phase for interactions with environment spins: φ = 0.0001; for details see [9,10].

The master equation with reset mechanism
Compared to (3), the master equation that models a gas-type system with reset mechanism has an additional term L reset , which we describe as follows. With some probability rδt particle i, i = 1...N , is reset during the time interval δt to some specific state |χ i . The other qubits are left in the state tr i ρ. The change in the density matrix during the time δt due to L reset is (δρ Observe that the time interval must be longer than the timescale of any of the involved processes but short enough so that we can replace it by the time differential to obtain, for the rate of changeρ = ∂ρ/∂t, the following master equation: Before we proceed to discuss the solution of (7) let us establish that the problem is well-defined, i.e., that the master equation leads to a completely positive map, which is true when the master equation is of Lindblad form. For the noise part this is known, so we have to bother only about the reset part. Since L reset = N j=1 r(|χ j χ j |tr j ρ − ρ) is local we have to show that each summand is of the form n ]) where the σ (j) s are Pauli operators and L j must be positive (semidefinite) matrices. We expand ρ and |χ j χ j | in the σ-basis as ρ = 3 k1,...,kN =0 a k1,...,kN σ q . We insert these expressions into L reset and also into the Lindblad-expression, collect the coefficients that belong to each σ-matrix of ρ using the scalar product, and compare the coefficients a kl in each expression, which leads us to a simple linear system of equations for the L j mn . Solving this system of linear equations we obtain The eigenvalues of L j are r/8, and eigenvalues r/8, 0, r/4. We note that also a mixed reset state would be fine to ensure that the L j are positive semidefinite. Because the sum of positive matrices is a positive matrix, and because we know that the noise terms also have positive L-matrices, we have shown that the master equation is of Lindblad form and preserves the positivity of the density matrix.
Up to this point, we have modeled interacting, gas-type systems coupled to a noisy environment. We have described a toy model consisting of only two particles by a master equation and compared predictions about the entanglement properties of steady states to simulations with a spin gas as an example for such gas-type systems. We have seen that in general there will be no entanglement in the steady state. We have extended the example of the spin gas by allowing particle exchange with a "cold bath" of particles in standard states (or an equivalent reset mechanism), and we have found that steady states of such systems can be entangled. We have derived a master equation that models systems with reset mechanism and have proved that the master equation is of Lindblad form. In the following subsection we study the solutions of (7).

Solution of the master equation for the gas-type model
In principle, the solution to the master equation (7) with noise channels as in equation (4)  Mapping ρ to a column vector C containing the 16 coefficients C 0000 , C 0001 , . . . , C 1111 of the density matrix and accordingly mapping the Liouville operator L to a 16 × 16matrix Λ we get the equivalent 16 coupled linear differential equationsĊ = ΛC with solution C(t) = e Λt C(0). To compute the matrix exponential we need the spectral representation of Λ, i.e., we must solve the eigenproblem ΛC λ = λC λ . Observe that the steady state (if it exists) is given by the eigenvector C 0 corresponding to the eigenvalue λ = 0.
In the following we will first analyze the solutions for an Ising interaction Hamiltonian and later generalize to generic cases.

Ising Hamiltonian.
We specialize to the Ising Hamiltonian (2) as (effective) interaction Hamiltonian, H = H Ising , and to a specific reset state, namely |+ +| for both qubits. The free Hamiltonian is as before. We can solve the problem through spectral decomposition of the Liouville operator L(ρ) = −i[H total , ρ] + L noise ρ + L reset ρ, but the expression for the corresponding matrix e Λt is very lengthy.
One obtains shorter expressions if one does not solve all 16 differential equations at once through the matrix exponential, but step by step, since not all differential equations are coupled. Still, we have chosen to move the solution derived in this way to Appendix A not to overburden the text with technical details.
For illustration, we will restrict the noise to the special case of a dephasing channel (5) in the following.

Dephasing channel
As pointed out above, the solution is given by the spectrum of the total Liouville super-operator defined byρ = Lρ and its corresponding eigenvectors. In the case of a dephasing channel one obtains the eigenvalues with multiplicities {1, 2, 1, 2, 1+1, 2+2, 2+2}, respectively. The eigenvector belonging to the eigenvalue 0 represents the density matrix in the steady state, and we will come back to this matrix in the next subsection. The full solution, derived by solving the differential equations in a step-by-step manner as explained above, is given in Appendix A.
To demonstrate the time-evolution of an initial density matrix governed by the master equation we plot the entanglement between the two qubits as a function of time. Note that Figure 3 is based on the analytic solution given in Appendix A, and the plot shows how the negativities of different initial states approach the final negativity of the steady state. We choose γ −1 as unit timescale (setting γ = 1). The parameters r, g, ω here have the fixed values 10γ, 5γ, 5γ, respectively. The initial states are weighted graph states [9,12] with density matrix U (ϕ)|+ +|U † (ϕ) where U (ϕ) = diag(1, 1, 1, e iϕ ). Through the parameter ϕ we can continuously tune the entanglement in the initial state from the product state |+ +| for ϕ = 0 to the maximally entangled, Bell-equivalent state for ϕ = π. States that are initially highly entangled are first driven into separable states before the steady-state entanglement value is approached from below. Vice versa, an initial product state gets highly entangled first, before the steady-state value for the entanglement is reached from above. As we said earlier, to display the full analytic solution of (7) for more general noise channels would be quite space-consuming. We will not present it since we are primarily interested in the entanglement properties of the steady state. In the following we will discuss these properties for the master equation with general local noise channels.

General steady states with Ising Hamiltonian.
As we have seen, any initial state of the density matrix evolves exponentially fast into a steady state on a characteristic timescale given by the largest non-zero characteristic exponent (or the smallest in absolute values since they are negative). The characteristic timescale thus depends on the parameters of the master equation, too. The steady state is also a function of these parameters.
To smooth the presentation, we have again transferred the steady-state solution of (7) with Ising Hamiltonian, local noise channels as in (4)   to Appendix B. Here, we illustrate the solution with a plot. In Figure 4 we choose B −1 as unit timescale and fix the values C = B, ω = 20B and the dimensionless parameter s = 0.1. Then, for certain values of r and g, measured in units of B, we see entanglement as measured by the negativity in the steady state. To get entangled steady states when g becomes small, we have to go to higher reset rates r. However, g cannot be arbitrarily small. There is a weak coupling threshold below which no reset rate can ensure entanglement in the steady state. This threshold depends on the parameters of the decoherence processes, and its existence is intuitively clear. If the decoherence processes simply dominate the entangling processes, then no entanglement can be created by any means. To see this better, we momentarily put ω = 0 for simplification. Let us also turn our attention once more to the dephasing channel as a special case of the noise terms of the quantum-optical master equation. Then, in the steady state, the anti-diagonal coefficients are all the same, namely , and the diagonal elements all have the values 1 4 . All other matrix elements are given by the Hermiticity of the density matrix. We compute from the above expressions for the density matrix the following analytic expression for the negativity in terms of the parameters g (Hamiltonian interaction), γ (strength of the dephasing channel), and r (reset rate): Equation (9) contains the full information about the entanglement properties of the two qubits. Note that N = N (g,r) depends, in fact, only on two parameters,g = g/γ andr = r/γ. In Fig. 5(a) we see a plot of the negativity function N. The key feature is the color-coded region in the r-g-plane with steady-state entanglement, where a darker color indicates higher entanglement. The entangled region is bounded by the red line given by one of the roots of the non-trivial part of equation (9). Outside of this region, the state is separable (white area). The entangled region approaches asymptotically the straights g = 2γ and g = r plotted black in Fig. 5(a). The asymptotic line g = 2γ is independent of r and simply tells us that, in the weak coupling regime, decoherence/noise will always triumph over the Hamiltonian part that tries to create entanglement as pointed out before. That is, as a necessary condition, we need to be above this threshold to observe entanglement. Three lines are marked in the colored, entangled region: (i) The upper, white line is the maximum in g-direction (at constant r).
(ii) The lower, white line is the maximum in r-direction (at constant g).
(iii) The middle line in black is the straight g = r/(1 + √ 3). To this middle line the upper and lower white curves go asymptotically for large g, r.
The global maximum of the negativity is on this middle line at infinity with a value of approximately 0.0915, about 20% of the maximally possible value. The darkest, most entangled area in our plot has negativity approximately 0.068. Fig. 5(b) shows a cut at g = 5γ through the color-plot. Most notable is the existence of a threshold value for r/γ above which entanglement is present in the steady state.

Steady state entanglement with and without reset.
In this part we want to compare steady-state entanglement that is due to the reset to entanglement that is due to a special combination of interaction Hamiltonian and decoherence process as in subsection 2.1. There, the interaction Hamiltonian is x while the free Hamiltonian are given by (1). The decoherence process (4) is determined by the parameter choice C = B/2. In 2.1, it has been established that the steady state is entangled, but only in a finite temperature range. If we add the reset mechanism to the master equation of this example we can show that the steady state is entangled for arbitrary temperatures. As reset states we choose the eigenstate |1 = −σ z |1 of the Pauli operator σ z .
(ii) There is a threshold reset rate r, above which an entangled steady state exists for arbitrary temperature.
From the coefficients of the steady-state density matrix (10) we also see that the entanglement created by the reset stems from a different density matrix than the entanglement present without reset. In Figure 6 this is visible in the region of small r, where the reset tends to destroy this latter entanglement. Then, for larger r, the effect of the reset mechanism kicks in.
x σ (2) x +ω/2(σ (1) z + σ (2) z ) as a function of the reset rate r and temperature-dependent parameter s. The noise is described by the quantum-optical master equation (4) with B = 2C = 1. With B −1 as unit timescale the other parameters are given by g = 2B and ω = 2B. The entanglement due to the reset mechanism exists for all temperatures, while the entanglement without reset mechanism (r = 0) vanishes above a certain temperature threshold. (b) Cut at constant s = 0 (orange curve, corresponding to zero temperature) and s = 0.5 (blue curve, corresponding to infinite temperature).
Up to now we studied rather special interaction Hamiltonians. In the following we demonstrate the genericity of entanglement that is present in a steady state due to a reset mechanism.

2.4.4.
Generic cases. In the previous example we have pointed out that entanglement, if present at all without reset, stems from a special combination of Hamiltonian and decoherence process. One may ask whether adding a reset mechanism with special reset states is not just as artificial as the choice of special combinations of Hamiltonians and decoherence processes. We are now going to show that one and the same reset with fixed reset states can lead to steady-state entanglement for many combinations of Hamiltonians and decoherence processes. The reasoning was the opposite without reset, where only very few combinations of Hamiltonians and decoherence processes lead to steady-state entanglement. We can also relax the condition that the reset states are pure states to a certain extent. Hence, we show that we have found generic features by generalizing the system in various directions.
(i) The qualitative behavior of the two-qubit model does not depend on the particular choice of the interaction Hamiltonian or details of the decoherence model other than its local action on individual qubits. Figure 7 shows e.g. steady-state entanglement for an XYZ Hamiltonian as function of reset rate r, and decoherence described by the noise operator L noise . The qualitative behavior is similar to Figure 5(b) or Figure 6, and we observe steady-state entanglement even for infinite temperature of the bath.
x + 0.3σ x )) + ω/2(σ (ii) The idealized reset mechanism we consider can be replaced by a more realistic imperfect reset mechanism. In this case, fresh particles are in mixed states with sufficiently low entropy rather than in pure states (with entropy 0). Still, the steady state turns out to be entangled. When we vary r there is a new, third threshold. First, for very small r, there can be an entangled steady state, which is not due to the reset and which is present only for a finite temperature range above zero. Second, there is one threshold value for r above which the steady state is entangled due to the reset mechanism (for arbitrary temperature). Third, whereas for pure reset states the entanglement goes down to zero again only in the limit r → ∞ (permanent projection to the product of the reset states), the entanglement goes down to zero for some finite r if the reset states are mixed states. This behavior is easy to understand, since it is more difficult for the interaction Hamiltonian to create entangled states from mixed reset states. The higher the entropies of the reset states are, the smaller is the range of the reset rate r for which there is entanglement in the steady state. This range can also become zero, so we must demand reset states of sufficiently low entropy. The picture that emerges from all these results is the following. Entanglement can prevail in dissipative, open quantum systems that are far away from thermodynamic equilibrium. For gas-type systems treated in this section, a reset mechanism can evoke steady-state entanglement even for infinite temperature of the environment generically, i.e., independently of the specific form of the interaction Hamiltonian or decoherence channel.
In the next section we show that steady-state entanglement appears also in strongly coupled systems with an appropriate reset mechanism.

Strongly coupled systems
In gas-type systems, we can treat the local noise channels separately for each qubit as explained above. If these local channels correspond to a heat bath, they drive each qubit individually to the thermal state of the local, free Hamiltonian, i.e., they populate the eigenstates of the free Hamiltonian according to the Boltzmann factor. Although the effective interaction Hamiltonian in the master equations is represented as continuously acting, the physical interaction process in gas-type scenarios is viewed as a short collision event. Hence, the interaction Hamiltonian does not influence the energy spectrum considerably and does not modify this thermal state. In strongly coupled systems, on the other hand, interactions of quanta of the heat bath with the system qubits affect the system as a whole. In this sense, the decoherence process acts globally on the system, inducing transitions between joint eigenstates. In this section we will shortly discuss the master equation describing a strongly coupled spin system in contact with a thermal, photonic bath. We will see that the resulting equilibrium state, the thermal state, can be entangled below a certain temperature threshold if the ground state of the Hamiltonian is entangled. When we add a reset mechanism we find that, in contrast to the gas-type scenario, entanglement in the steady state can exist only below a certain temperature threshold. However, the novel feature is that this threshold is typically much higher than for the thermal state. Finally, we describe the influence of the master equation parameters on the respective steady state, and, with this insight, formulate a general condition under which a reset mechanism can lead to an entangled steady state.
Let |a (|b ) be momentary eigenstates of some non-degenerate system Hamiltonian H(t) with eigenenergies ω a (ω b ) ‡. We define N ab := (e β(ωa−ω b ) − 1) −1 with β = 1/T being the inverse temperature. Often N ab is written asn, and we explained the connection to the parameter s in the last section. The master equation for a spin system, coupled with strength γ to a heat bath consisting of photons, is [13] Here, g(ω) is the spectral density, for which g ba = g(ω b − ω a ) if ω b > ω a and g ba = 0 else. For small system one can justify to treat the spectral density as constant (g = 1 if ω b > ω a ) and merely tune the overall coupling constant γ [13]. Observe that we did not include non-radiative contributions as opposed to the master equation (3) with noise terms (4). The master equation (12) drives any initial density matrix to the thermal state of inverse temperature β. That means, the ground state and also the excited states are populated according to the canonical distribution.
We will study the master equation (12) for an Ising Hamiltonian with transverse magnetic field, briefly discuss the solution without reset mechanism (thermal state), and then turn to an analysis of the full master equation with reset mechanism.
where N = (2 + 1/2|(−1 + √ 1 + 4b 2 )/b| 2 ) −1/2 provides normalization. We exclude the case b = 0 where the ground state would be degenerate, a case not properly described by the master equation (12). The ground state of this system is the first eigenvector in equation (14). Since this state is entangled, so is the thermal state below a certain temperature threshold. We see this directly from the negativity of the thermal state ρ thermal = exp{−βH}/ tr exp{−βH}, which is For any fixed g and b, the non-trivial part in this formula goes to the value −1/4 when β → 0, while it goes to 1/(2 √ 4b 2 + 1) > 0 for β → ∞. The threshold value for β, where the non-trivial part becomes exactly zero, can be easily computed numerically for any given parameters g and b. In terms of β, the thermal mixture of the eigenstates is separable below this threshold, which, in terms of T = β −1 , means that the mixture is separable above that critical temperature. From 15 one can see that the critical temperature grows linearly with g and monotonously, but sub-linearly with b.

Master equation with reset
We keep the Ising Hamiltonian with transverse magnetic field as above, but extend the master equation (12) by the reset term with N = 2 and reset state |χ i = |+ for both qubits. We solve the resulting master equation numerically. In Figure 9 we see how the entanglement for different reset rates r develops over time t from the value zero in the initial product state ρ = | + + + + | to its final value in the steady state while all other parameters are kept fixed (see figure caption). We notice that small, non-zero reset rates decrease the entanglement in the steady state until it is gone, while larger rates can bring entanglement back.

Influence of the parameter r.
To explain this effect, recall that the master equation mimics the averaged density matrices that would be obtained from (infinitely many) simulation runs of the system. The reset rates of the master equation are related to probabilities that in a simulation a reset took place during a certain time interval. Although our reset processes are strictly speaking local, let us assume for the sake of argumentation that the reset happens on both qubits simultaneously, thus effectively restarting the process again from the beginning whenever a reset occurs in a simulation. For small rates r, i.e., for small probabilities that a reset takes place in the simulation, the system can come close to its thermal equilibrium state before it is reset. When we average the density matrices over many simulation runs, we average matrices that are mostly close to the unique thermal equilibrium state, and hence also the mixture will still retain entanglement. When the rates get larger, the density matrices over which we average become more and more diverse since they will be far from equilibrium and fluctuations occur. As a consequence the average density matrix will have no entanglement. When the reset rate is above a certain threshold, we will find entanglement in the system again (as we did in the gas-type systems) because now the density matrices over which we average become similar again. Now, they are close to the state that has unitarily evolved for a time of order 1/r from the initial reset state. In the limit r → ∞ the state is constantly kept in the initial product state with zero entanglement. In this way we can understand how the two entangled regions arise. The first is an artifact of the entangled thermal state that is more and more destroyed by the reset mechanism. This first region could also be present in a gas-type model. The second region is the one that is really created by the reset mechanism just as in the gas-type model. We can directly see in Figure 9 that entanglement in a thermal state is a truly different effect from entanglement that is created by the reset mechanism. Figure 9. (a) Solution of the master equation (12) with reset term (16) for a strongly coupled system. The unit timescale is γ −1 and the time t is plotted in these units. The temperature parameter was chosen as β = 1000, and, at this low temperature, the steady state of (12) is entangled even without reset, the other parameters being g = 10γ and b = 0.1. When increasing the parameter r of the reset mechanism, the steady-state density matrix changes, as explained in the text, fist becoming separable and then entangled again. (b) Cut through the same plot for different r.

3.2.2.
Influence of the parameter γ. Although the two effects are truly different, this does not mean that for certain parameter regimes, the two regions cannot overlap, see Figure 10. Imagine the coupling to the photon bath, γ, is increased. This means that the system is driven towards thermal equilibrium faster than before. Hence, following the arguments from above, the system can tolerate higher reset rates before the entanglement in the first region is destroyed. The stronger coupling to the photon bath suppresses the entanglement in the second region, and as an overall effect we see that the two regions need not be separate. Note that the entanglement in the thermal state for r = 0 is independent of γ because then it does not matter how fast equilibrium was approached.  (and will approach the product state | − − for b → ∞). Hence, at zero temperature, the entanglement in the equilibrium state for r = 0 will go down for increasing b (see red line in Figure 11(a)). For a thermal state with T = 0, i.e., β = ∞, the ground and exited states get mixed. When b is small, the splitting between ground and first exited state is small, and the mixture will be close to a separable state. For increasing b, the larger energy split leads to an increased population of the ground state relative to the exited states at the same temperature 1/β, and the entanglement will increase. When b gets even larger, the thermal state will be close to the ground state, but we know, that the ground state for large b is only weakly entangled. Hence, there will be some b for which the entanglement in the thermal state is maximal (see red line in Figure 11(b), and the identical line in Figure 11(c)). When we switch on the reset mechanism, we see that an increasing r destroys the entanglement in the thermal state as before, but for increasing b the regions one and two (artifact of thermal state vs. true reset state entanglement) quickly overlap. That is because the magnetic field tends to drive the state towards | − − , whereas an increasing reset mechanism drives the state towards | + + . If we choose the reset state as |− , the two effects do not compete and both drive the state towards an unentangled product state (see light grey areas in Figure 11(c) and compare to the light grey areas in Figure 11(b)). How fast an increasing r destroys the entanglement in the thermal state is almost independent of the entanglement in the thermal state. For larger r the influence of b plays less and less a role and there is almost no difference between Figures 11(b),(c) for large r.

Influence of the parameters g and β.
The coupling strength of the Hamiltonian g also splits the energy levels. Hence with increasing g there will be more entanglement in the thermal state at some finite temperature. Again, the speed with which an increasing reset rate destroys the entanglement in the thermal state is almost independent of g (see Figure 12(a)). Most interesting is the influence of the temperature 1/β. Figures 12(a)-(c) show plots which contain information about the equilibrium-state entanglement for different temperatures (r = 0). We see how the thermal states get less and less entangled for increasing temperatures, so that region Figure 11. Influence of the relative magnetic field b. The unit timescale is given by γ −1 and g = 50γ. Plots (a) and (b) show the difference between (almost) zero temperature (β = 10000) and some finite temperature (β = 0.2). The behavior of the red curves, for r = 0, are explained in the text. Plots (b) and (c) demonstrate the influence of different reset states in connection with the magnetic field b. Plot (b) has reset state |+ , (c) has reset state |− . one vanishes quickly as expected. But, for the same temperatures, the reset rate r can still produce entanglement in a steady state! On the other hand, for fixed g, there is some temperature threshold above which no reset rate can produce entanglement. This threshold is in contrast to the case of gas-type systems. Figure 12. Negativities for increasing temperature (decreasing β = 1/T ). The unit timescale is given by γ −1 , the relative magnetic field by b = 0.1. The parameter region in the r − g plane, where steady-state entanglement occurs, becomes smaller for increasing temperatures. However, the temperatures, for which entanglement can exist with reset mechanism, are much higher than the temperatures, for which the thermal state is entangled. The red line in (c) corresponds to the red line in Figure 13 (see that figure caption and the text) and is drawn for comparison only.

General conditions for steady-state entanglement
This discrepancy between gas-type and strongly interacting systems raises the deeper question: What are the conditions under which the reset mechanism can create entanglement in the steady state? We already see that the reset mechanism is different from a cold bath -equivalently the replacing of system particles by fresh, standard ones -because we would expect that some cold bath can always counteract the influence of the hot photon bath. The condition that the solution of the master equation (a completely positive (CP) map) at r = 0 is entangling at some point in time is certainly necessary, since, as pointed out before, the reset mechanism does not introduce entanglement itself. The question is: Is a solution of the master equation at r = 0 that creates entanglement on some short time scale sufficient such that the solution for some r > 0 is a CP-map with entangled steady state? Unfortunately this is not true. The reset rate r itself can influence the solution of the master equation in such a way that although the solution for r = 0 was entangling on some short time scale, the solution for larger r need not be.
The condition for steady-state entanglement is the following. If the solution of the master equation for some r > 0 is entangling on a time scale of order 1/r then the steady state of this solution will also be entangled. For illustration, we think once more in terms of many hypothetical simulation runs. As stated above, the states over which we have to average will be close to the state that has unitarily evolved for a time 1/r from the initial reset state. When this state has a certain amount of entanglement, then so does the mixture of states close to it (see Figure 13). Figure 13. Illustration of the condition under which a reset mechanism can create entanglement in the steady state. At time γt = 2 the state is already close to the true steady state, so the red line corresponds to the red line in Figure 12(c) since g = 10γ in this plot, and the other parameters are the same as in Figure 12(c). The reset can create entanglement in the steady state if, for given r, there is entanglement in the time-evolved state at time t ∝ 1/r (see text for details). The green curve is the curve t = 2/r illustrating this result.
The condition does not only hold for the specific Hamiltonian or the specific heat bath chosen here. It is valid for a large class of Hamiltonians (with appropriately chosen reset states) and baths.
In the next section we treat the multipartite case.

Multipartite case
For multipartite spin systems, we will show that the reset mechanism can create steady-state entanglement in a similar way. The parameter regions where this happens are comparable to the 2-qubit case. The values of the negativity, or rather its generalization, the average negativity, stay almost constant with increasing system size. Note, however, that larger systems could have larger negativities, so if we divide the actual (constant) negativity by the maximal possible negativity, then this quantity would go down for growing system size. We will also consider entanglement in reduced density matrices. Since in a reduction from N to, say, 2 qubits the traced out N − 2 qubits act as an additional noise source, it is not surprising that the parameter region where we find steady-state entanglement in the reduced systems shrinks with growing N . If the number of particles fluctuates according to some distribution, our best description of a reduced density matrix is a mixture of reductions originating from different system sizes. It is remarkable that even in this case there is some parameter region, where steady-state entanglement is found.
In the following, we motivate and explain the entanglement measures we are using and then demonstrate the above features in both gas-type systems and strongly interacting systems.

Gas-type systems
It is straightforward to generalize the master equation for the gas-type system, equation (7) to N > 2 qubits. The relevance of the entanglement quantities we are going to use needs to be motivated, though. We turn once more to the example of a spin gas. Imagine that the spin gas is in a box of volume V . This box has one semi-permeable wall through which particles can leave the box, and another through which particles from a "cold reservoir" can enter. By "cold reservoir" we mean that the quantum state of the particles in this reservoir is in some sufficiently pure standard state. The motional degrees of freedom of theses particles, on the other hand, are in thermal equilibrium with the outside environment just like the system particles in the box. Assume that the density of the gas in the box is η. Then there are on average ηV =: λ particles in the box. The distribution of the number of particles that are in the box is a Poissonian p λ (n) = e −λ λ n /n!. When we observe the spin gas after certain time intervals, which should be long enough such that the gas always reaches its equilibrium state, we sample the distribution and get information about the density matrices with a corresponding number n of qubits. The density matrices we can reconstruct after we collected a certain, sufficient amount of information is close to the steady-state density matrix of the master equation for n qubits. We are interested in the entanglement properties of the gas and we will look at different aspects of entanglement in the following.
All these aspects of entanglement are quantified by measures that are based on the negativity or the average negativity. The average negativityN is the negativity averaged over all possible bipartitions of the system [9]. Non-zero average negativity ensures the presence of some form of entanglement in the system.
Specifically, we study three types of entanglement: (i) The average negativity of the density matrices with n qubits, averaged over the Poissonian distribution of the number of particles in the system: N (ρ n ) p λ (n) . (ii) The negativity of reduced 2-qubit density matrices averaged over a renormalized, truncated Poissonian distribution: N (ρ n→2 ) p λ (n≥2) .
Now, we lay out, what these quantities mean, and which aspect of entanglement they describe.
(i) If we ask how much entanglement we find in the system on average we are led to the quantity N (ρ n ) p λ (n) , which is the expectation value of the average negativity of density matrices with different n, where p λ (n) is the Poissonian probability distribution for the number n of particles in the system. Observe that if we disallowed the fluctuation of the particle number in the system and introduced the reset mechanism by other means (e.g. measurement and decay to standard state inside the box), the quantity of interest would simply beN (ρ N ) for fixed system size N .
(ii) When we look at subsystems we are led to slightly different quantities. Let us fix the subsystem size to 2 qubits. We call the reduced density matrices of originally n qubits ρ n→2 where we assume n ≥ 2. In gas-type systems there can also be zero or one particle in the box (especially if η ∝ λ is small), and the entanglement is simply zero in these cases. Since a "reduction" of a one or zero-qubit density matrix to a 2-qubit density matrix makes no sense, we simply exclude these cases and rescale the truncated Poissonian p λ (n ≥ 2) to the distributionp λ (n ≥ 2). The quantity N (ρ n→2 ) p λ (n≥2) therefore tells us how much entanglement a subsystem of two qubits contains on average (for 2 qubits N =N ).
(iii) When we look only at a subsystem of two qubits disregarding the number of particles n in the system, then our best description of the 2-qubit density matrix is the average density matrixρ := ρ n→2 p λ (n≥2) with entanglement N (ρ).
When we compare the three kinds of entanglement defined above we see that the conditions that one of these quantities be non-zero are increasingly stringent. To find entanglement in the reduced system is a more stringent condition since tracing out the other particles has the same effect as an additional noise source. Also, to find entanglement in the averaged density matrixρ is a stricter condition since the averaging increases entropy, i.e., tends to make the matrix more mixed. If we keep the particle number fixed, N ≥ 2, there is just the quantity N (ρ N →2 ) to describe the entanglement in the reduced state.
To compute the average negativity is a hard task. The system size, and hence the number of differential equations we must solve, and also the number of bipartitions scale exponentially. To simplify the computation, we consider a symmetric situation, where all qubits interact pairwise via Ising interactions and are subject to dephasing noise (5). In Figure 14 we plot the three entanglement measures N (ρ n ) p λ (n) (blue, dashed curve), N (ρ n→2 ) p λ (n≥2) (orange, solid curve), and N (ρ) (black, dasheddotted curve) for λ = 2 and for a qubit number that fluctuates between two and five. Since the meaning of these measures is different, one cannot compare the absolute values represented by the curves directly with one exception. The points, where the curves become non-zero, must, from left to right, appear in the order explained in the previous paragraph, i.e., blue first, representing measure (i), orange next, representing (ii), and black last, representing (iii).  (1) is ω = 50γ. Kinks in the orange curve stem from the averaging over negativities whith different supports (see (ii) in the text).

Strongly coupled systems
Eventually, we study equations (12) and (16) in the multipartite case. We obtain similar results as in the case of gas-type systems underlining again how generic the reset mechanism is. Although a fluctuation of particles may seem less natural in strongly coupled systems as compared to gas-type systems, we will look at the exact same entanglement quantities, so that the corresponding plots are directly related to each other. As Hamiltonian, we choose a sum of pairwise Ising interactions and magnetic fields in x and z direction according to where the small gradient magnetic field in z direction is introduced for technical reasons to lift degeneracies in the Hamiltonian. Figure 15 shows a plot of the same entanglement measures as in the previous subsection. As in the gas-type scenario, the feature that entanglement can be created by a reset mechanism holds also in the multipartite case. While we have shown the genericity of the reset mechanism with respect to the Hamiltonian and noise process already in previous sections, here we demonstrate that the reset mechanism is also generic with respect to system size.

Summary
We mechanism that "resets" the particles, at a certain rate, into a single-particle, lowentropy state. For a 2-qubit toy model of a gas-type system, we have analytically solved the master equation consisting of a Hamiltonian part, a noise channel, and the proposed reset mechanism. For special cases we have been able to give closed expressions for the entanglement as a function of the parameters of the master equation. We have extended the analysis to similar models with other interaction Hamiltonians, decoherence models, and imperfect reset mechanisms. We have treated the situation of strongly correlated systems by the same means and we have given conditions under which steady-state entanglement arises in this case. Finally, we have shown that in systems consisting of more qubits, and even in systems with fluctuating particle number, steady-state entanglement can prevail.
Many systems are conceivable for an experimental realization of such dissipative, open quantum systems. For instance, one may consider ions in microtraps that interact via an induced dipole moment [14] leading effectively to a continuously operating Ising interaction. Decoherence, dominated by dephasing noise, appears naturally in such systems, and the reset mechanism may, e.g., be achieved by periodically applying a πpulse that couples the internal level |1 to a metastable auxiliary level |a that decays rapidly to |0 . The state afterwards is always |0 , which can be mapped to |+ by a subsequent Hadamard operation.
For charge manipulated quantum dots, the exchange interaction leads to a continuously operating Heisenberg interaction between neighboring electron spins by lowering the potential barrier [15]. The effect of surrounding nuclear spins may be described by dephasing noise, while the reset mechanism can consist in replacing an electron by a fresh one from the surrounding Fermi sea, prepared in a suitable state (e.g., |0 ).
Atomic beams interacting via a cavity mode [16] may also serve as a toy example of such systems far from thermodynamic equilibrium.
Note that a reset mechanisms could be realized in many physical ways, including a measurement with subsequent preparation of the state, coupling or decay to metastable auxiliary states, as well as replacing a qubit by a fresh one. While these suggestions for implementations aim at demonstrating how a reset mechanism would be realized in experiments, the effect itself is generic and other realizations are conceivable. In particular, one might try to find such a reset mechanism in less controlled, maybe even biomolecular systems consisting of many particles. themselves and to the diagonal, so we can solve them next. The same is true for the pair C 0010 and C 0111 . Finally, the anti-diagonal coefficients C 0011 and C 0110 are coupled to C 0001 ,C 1011 ,C 0010 , and C 0111 , or to their complex conjugates. We solve these as a last step, and all other coefficients are given by the Hermiticity of the density matrix. The solution is now straightforward in principle. However, the expressions for the matrix coefficients are still space-consuming, so we will give them only for the special case of the dephasing channel.
Furthermore, C 0010 = C 0001 , and C 0111 = C 1011 are obtained from C 0001 by replacing g → −g, s → (1 − s) in the numerator. The other coefficients are given by the Hermiticity of the density matrix. Observe that the dephasing and the depolarizing channels are included as special instances of the parameters B, C, s in this analytic expression. The dephasing channel results from putting B = 0, s drops out, and renaming C = 2γ. The depolarizing channel is given by s = 1/2, B = C, and renaming C = 4γ/3. The plot in Figure 4 is based on this solution.