Quantum master equations for entangled qubit environments

# Quantum master equations for entangled qubit environments

Shakib Daryanoosh Centre for Engineered Quantum Systems, Department of Physics and Astronomy,
Macquarie University, Sydney, NSW 2122, Australia
Ben Q. Baragiola Centre for Quantum Computation and Communication Technology, School of Science,
RMIT University, Melbourne, Victoria 3001, Australia
Centre for Engineered Quantum Systems, Department of Physics and Astronomy,
Macquarie University, Sydney, NSW 2122, Australia
Thomas Guff Centre for Engineered Quantum Systems, Department of Physics and Astronomy,
Macquarie University, Sydney, NSW 2122, Australia
Alexei Gilchrist Centre for Engineered Quantum Systems, Department of Physics and Astronomy,
Macquarie University, Sydney, NSW 2122, Australia
July 16, 2019
###### Abstract

We study the Markovian dynamics of a collection of quantum systems coupled to an irreversible environmental channel consisting of a stream of entangled qubits. Within the framework of repeated quantum interactions, we derive the master equation that describes the dynamics of the composite quantum system. We investigate the evolution of the joint system for two-qubit environments and find that (1) the presence of antidiagonal coherences (in the local basis) in the environment is a necessary condition for entangling two remote systems, and (2) that maximally entangled two-qubit baths are an exceptional point without a unique steady state. For the general case of -qubit environments we show that coherences in maximally entangled baths (when expressed in the local energy basis), do not affect the system evolution in the weak coupling regime.

###### pacs:
03.65.Yz, 03.65.Aa, 42.50.Dv, 42.50.-p,03.65.Ud

## I Introduction

Open quantum systems are the subject of extensive research since physical quantum systems cannot be entirely isolated from their surroundings. The influence of the environment often manifests as unwanted noise that can thwart attempts to exploit intrinsic quantum properties for quantum computing, communication, and metrology Barnett et al. (2017). Open systems tend to lose their key quantum properties—coherence and entanglement—as they interact with the environment. This is not inevitable, however, and much research has focused on engineered environments Müller et al. (2012) for various tasks including quantum computing Verstraete et al. (2009); Pastawski et al. (2011); Leghtas et al. (2015) and the generation of novel steady states Lin et al. (2013); Kienzler et al. (2014); Baragiola and Twamley (2018). This can be achieved through a combination of precise structuring of system-bath interactions and preparation of the environment in particular states.

One can take many approaches to the description of open quantum systems. In standard quantum optical treatments, the electromagnetic field serves as the environment Carmichael (2008); Wiseman and Milburn (2010); Breuer and Petruccione (2003), and the dynamics of the reduced quantum states is given by a master equation (ME) after the environmental degrees of freedom are traced out. An alternative approach is that of repeated quantum interactions Brun (2002); Pellegrini (2010); Bruneau et al. (2014), also called collision models Rau (1963); Scarani et al. (2002); Giovannetti and Palma (2012); Lorenzo et al. (2017); Francesco (2017), which treats the system-environment coupling discretely. The environment is comprised of a chain of identical and independent quantum ancillae which sequentially couple to the system and are then traced out. Taking a continuous limit Attal and Joye (2007), the resulting dynamical map on the reduced system state becomes a ME. This formalism has been applied in the context of quantum thermodynamics Horowitz (2012); Lorenzo et al. (2015); Daǧ et al. (2016); Strasberg et al. (2017); Li et al. (2018), as well as quantum optics and information Gross et al. (2017). The framework of repeated quantum interactions has also proven useful for the study of correlated quantum channels Korotkov (2002); Giovannetti and Palma (2012); Rybár et al. (2012). When the ancillae are entangled, the reduced-state dynamics can exhibit non-Markovian behavior Kretschmer et al. (2016); Lorenzo et al. (2017), such as propagating single-excitation states Dabrowska et al. (2017).

In this article we consider a correlated environment that interacts with many quantum systems simultaneously. The environment consists of a stream of ancilla qubits, each coupled to its own system. The qubits arrive entangled with one another at each interaction time, but they are not entangled across different times, which allows us to derive a Markovian ME for the joint state of the quantum systems. Depending on the state of the qubits, the ME can generate non-trivial dynamics of the reduced systems. This provides the tools to tackle the problem of converting coherences and/or entanglement in the environment into quantum correlations in the system Wu et al. (2017); Černoch et al. (2018); Riedinger et al. (2018). Accordingly, we analyze in detail the pedagogical case of qubits in the bath, which can be prepared as a stream of entangled states, as it provides the canonical method for transferring qubit entanglement to system entanglement. When the bath is prepared arbitrarily close to a maximally entangled state, the system is driven to an entangled pure state, such as a two-mode squeezed state for the case of two optical cavities. Surprisingly, if the qubit bath is exactly in a Bell state, the system fails to converge to a unique steady state. For -qubit baths , we find that for certain multi-qubit entangled environments such as baths prepared in -states (which have nonzero elements only on the diagonal and antidiagonal entries of the state matrix when expressed in the local energy basis) Weinstein (2010); Hashemi Rafsanjani et al. (2012), the bath entanglement cannot be transferred to the systems, as a direct consequence of the weak-coupling limit

The manuscript is organized as follows. In Sec. II the underlying framework behind our analysis is explained. We first present and interpret the two-qubit bath ME in Sec. III. We give two forms for the ME, which are useful for pure-state and general mixed-state baths, respectively. In Sec. V and Sec. IV we apply the formalism to study the dynamics and steady states for two sets of quantum system: optical cavities and two-level atoms. The general case of an -qubit bath is dealt with in Sec. VI. Finally, in Sec. VII we summarize the findings of this paper and propose future directions.

## Ii Repeated interaction with a bath of n entangled qubits

We derive dynamical maps and master equations within the structure of repeated quantum interactions Brun (2002); Bruneau et al. (2014); Strasberg et al. (2017); Gross et al. (2017). In this formalism, a quantum system in Hilbert space couples to an environment which is comprised of a stream of identical and independent quantum systems such that . We assume the environment has infinitely many elements, although in principle it can be finite. Each environmental element sequentially interacts with the system over a short time interval of duration after which it no longer interacts with the system. Tracing over the environmental degrees of freedom yields a map on the system of interest in . This is similar to the standard scenario for open systems in quantum optics where a bosonic probe field interacts with the system in a continuous-in-time manner Wiseman and Milburn (2010). However, the situation here is different in two ways. First, the system-environment coupling is fundamentally discrete, although we will ultimately consider short-time interactions and take a continuous limit Attal and Pautrat (2006); Attal and Joye (2007); Pellegrini (2008). Second, the environmental systems are qubits rather than bosonic modes. This serves not only to model physical situations where streams of qubits interact with a fixed quantum system Nogues et al. (1999), but in addition the results fit into the framework of quantum computing and simulation Gross et al. (2017).

For each interaction time interval the total Hilbert space of the system plus the segment of the environment interacting at that time is , and the Hamiltonian is

 ^H=^H0+^H(l)I, (1)

corresponding to the bare Hamiltonians of the system and of the environment, , and the the system-environment interaction, . We take the environment to be comprised of a stream of entangled qubits with each qubit coupled to its own quantum system, see Fig. 1. That is, in each time interval , entangled qubits interact with quantum systems, which themselves are left arbitrary and can in general be remote from each other. Note that the bath qubits within a single time interval are entangled, but they are not entangled between time intervals. This type of environmental entanglement drives non-Markovian dynamics Korotkov (2002); Gough et al. (2011); Baragiola et al. (2012); Dabrowska et al. (2017); Baragiola and Combes (2017) and will be treated separately.

Each subsystem interacts with its respective qubit via a coupling operator . The bare and interaction Hamiltonians in the rotating wave approximation are

 ^H0 = n∑ℓ=1(ωℓ^c†ℓ^cℓ+ω\tiny Eℓ^σ†ℓ^σℓ), (2) ^H(l)I = n∑ℓ=1λℓ(^cℓ^σ†ℓ+^c†ℓ^σℓ), (3)

where and are the respective transition frequencies of the subsystems and the bath qubits, is the coupling strength between the -th subsystem and its qubit, and the bath qubit lowering operator is . In the interaction picture with respect to the bare Hamiltonian, , the joint unitary evolution is generated by the time-dependent Hamiltonian,

 ^H(l)I(Δt)=n∑ℓ=1λℓ(^c†ℓ^σℓe−iδℓΔt+H.c.), (4)

where is the detuning. The detuning is included here for completeness; henceforth, we focus our attention to resonant system-qubit interactions, . From the interaction, one models a specific reservoir by selecting a particular state for the environmental qubits. Investigating situations where the bath qubits are entangled is the focus of this article.

A dynamical map for the joint state of the fixed quantum systems is found by tracing out the the environmental qubits after the interaction, , generated by the Hamiltonian in Eq. (4). At each time interval the incoming bath qubits and the quantum systems are assumed to be in a product state, so the dynamical map is completely-positive and trace-preserving. Assuming that only one qubit interacts with each system in an interaction time and the joint state of the -qubit environment interacting in the time interval is , the dynamical map is given by

 ^ρ(tl+1)=TrE[^U(l)I(^ρ(tl)⊗^ρ(l)E)^U(l)I†]. (5)

The fact that each set of qubits, described by , is independent from other sets means that this dynamical map is Markovian (i.e. arises from a memoryless environment). This is because the joint state of the input environment is a tensor product state across interaction time intervals, Baragiola and Combes (2017). Thus, the dynamical map at every subsequent time interval is of the form of Eq. (5), using the system state from the previous time interval and a fresh environmental state . From now on, for the ease of notation we drop the explicit superscripts for the system and environment state unless confusion could arise.

## Iii Master equations for two-qubit baths

In this section we focus on two-qubit baths as the quintessential extension of the single-qubit baths that are typically studied Lorenzo et al. (2015); Strasberg et al. (2017); Gross et al. (2017); Dabrowska et al. (2017); Li et al. (2018). Two-qubit baths can exhibit nonclassical correlations including maximal entanglement. We investigate how two-qubit baths modify correlations between two remote subsystems—optical cavities in Sec. IV and two-level atoms in Sec. V. Rather than using the discrete dynamical maps in Eq. (5), we take a continuous-time limit and describe the reduced-state dynamics by a Markovian ME.

A ME can be derived from Eq. (5) under a set of standard assumptions. First, each environmental qubit spends the same amount of time interacting with its local system. Second, the rotating wave approximation, which has already been made in Eq. (4), requires that the interaction time is long in comparison to the system’s characteristic time . Third, the Markov condition requires that the -qubit environment in each time interval is independent of other intervals. Finally, the system-bath coupling is weak while . In the weak-coupling regime the unitary time-evolution operator is expanded up to second order in ,

 ^U(l)I = e−i^H(l)IΔt (6a) = ^1SE−i^H(l)IΔt−12(^H(l)I)2Δt2+O(Δt3). (6b)

Let us assume that in each time interval the two bath qubits are prepared in the pure state , where

 |ψE⟩=bee|ee⟩+beg|eg⟩+bge|ge⟩+bgg|gg⟩, (7)

and the coefficients satisfy

 (8)

We insert Eq. (6b) into the dynamical map for the reduced system state, Eq. (5), and then evaluate the terms under the two-qubit bath state in Eq. (7). The ME in Eq. (9), which describes the continuous limit of repeated quantum interactions, then arises in the limit of infinitesimal interaction, . This calculation, whose details can be found in Appendix A, yields the following master equation (, for the reduced state :

 ˙^ρ(t) =−i[^Heff1+^Heff2,^ρ]+4∑m=1D[^Lm]^ρ (9)

where the effective Hamiltonians are

 ^Heff1 = λ1(bggb∗eg+bgeb∗ee)^c1+H.c., (10a) ^Heff2 = λ2(begb∗ee+bggb∗ge)^c2+H.c., (10b)

the jump operators are

 ^L1= √γ1bgg^c1+√γ2bee^c†2, (11a) ^L2= √γ1bee^c†1+√γ2bgg^c2, (11b) ^L3= √γ1bge^c1+√γ2beg^c2, (11c) ^L4= √γ1beg^c†1+√γ2bge^c†2, (11d)

and the Lindblad superoperator is defined as

 D[^o]^ρ=^o^ρ^o†−12{^o†^o,^ρ}+, (12)

with . The relative rates are given by .

The master equation generates both coherent dynamics and incoherent dynamics in the reduced system state . The Hamiltonian terms in Eq. (10) arise from coherences within each separate bath qubit. The coherent, unitary dynamics they generate are analogous to coherent driving Wiseman and Milburn (2010). Simultaneously, the quantum system undergoes correlated dissipation as described by the jump operators in Eq. (11). Each jump operator drives a dissipative process given by combinations of loss () and heating () across subsystems 1 and 2.

Interestingly, the jump operators (11) are determined by the state in two two-dimensional subspaces of the qubit bath, spanned by either or . Equivalently, each subspace is spanned by two Bell states. Populations and coherences within subspace contribute to jump operators and , and similarly populations and coherences within subspace contribute to jump operators and . Coupling between the subspaces is due to certain single-qubit coherences that manifest in the effective Hamiltonians in Eq. (10).

### iii.1 Unentangled bath qubits

An elementary situation is the case of unentangled bath qubits, . Consider uncorrelated bath qubits, both in the ground state, . The result is the typical ME for independent loss on each subsystem through the jump operators and , with . The situation becomes more complicated if either bath qubit contains excitation. In this case the ME can involve all four jump operators in Eq. (11) even when the two bath qubits are uncorrelated. Consider the case where each qubit is in the state with , and with no initial subsystem correlations, . The ME becomes

 ˙^ρ(t)= −i[^Heff1,^ρ1]⊗^ρ2−i^ρ1⊗[^Heff2,^ρ2] (13) −√γ1γ2[^Heff1,^ρ1]⊗[^Heff2,^ρ2],

with effective Hamiltonians (for ),

 ^Heffℓ = λℓαℓβ∗ℓ^cℓ+H.c. (14)

The first two lines describes the independent evolution of each subsystem, as would be expected if one were to separately derive and add together the ME for each subsystem. Here the ME contains an additional term in the final line. It describes a classically correlated coherent driving of the two subsystems, but it does not generate entanglement and preserves the separable structure of the joint state, . Further, tracing over either subsystem gives the expected single-subsystem ME. Note that this term does not appear in typical bosonic MEs due to vanishing bath correlation functions.

### iii.2 Bath qubits in a general state

While the master equation in Eq. (7) was derived for a two-qubit bath in a pure state, it is also valid for the environment in a general two-qubit state,

 ^ρE=∑j,k∑j′,k′bjk,j′k′|jk⟩⟨j′k′|, (15)

with indices . In order to show the correspondence we transform the ME into another useful form by expanding Eq. (9) and collecting terms according to coefficients, we get

 ˙^ρ(t) =−i[^Heff1+^Heff2,^ρ] (16) +2∑ℓ=1γ↓,ℓD[^cℓ]^ρ+2∑ℓ=1γ↑,ℓD[^c†ℓ]^ρ +γ↓↓S[^c1,^c2]^ρ+γ∗↓↓S[^c†1,^c†2]^ρ +γ↓↑S[^c1,^c†2]^ρ+γ∗↓↑S[^c†1,^c2]^ρ.

We have defined a superoperator,

 S[^o1,^o2]^ρ\coloneqq^o1^ρ^o2+^o2^ρ^o1−12{^o1^o2+^o2^o1,^ρ}+, (17)

that is symmetric in the arguments. The coefficients in Eq. (16) are given by

 γ↓,1 = γ1(∣∣bgg∣∣2+∣∣bge∣∣2), (18a) γ↑,1 = γ1(|bee|2+∣∣beg∣∣2), (18b) γ↓,2 = γ2(∣∣bgg∣∣2+∣∣beg∣∣2), (18c) γ↑,2 = γ2(|bee|2+∣∣bge∣∣2), (18d) γ↓↓ = √γ1γ2bggb∗ee, (18e) γ↓↑ = √γ1γ2bgeb∗eg. (18f)

The first four coefficients in Eqs. (18a-18d) are positive and can be interpreted as rates. The final two Eqs. (18e) and (18f) are complex and the superoperator terms that they multiply may in general interfere with other terms including the local dissipators. This is indeed a consequence of Eq. (16) not being in diagonal form with respect to the jump operators. For a general state of the two-qubit bath, Eq. (15), we replace the coefficients in Eq. (18) according to

 bjlb∗j′l′→bjl,j′l′,

and use these coefficients in the ME given in Eq. (16). The MEs in Eq. (9) and Eq. (16) are identical, but each may be more useful for certain calculations.

## Iv Two remote cavities: two-mode squeezing

A stream of qubits interacting with a harmonic oscillator is the prototype for a variety of tasks. For example, two-level Rydberg atoms interacting with an ultra-high-finesse microwave cavity, have been used for quantum nondemolition measurements of photon number Nogues et al. (1999) and stabilization of Fock states in the cavity Sayrin et al. (2011). We consider here the natural extension of this system to two remote cavities Campagne-Ibarcq et al. (2018); Didier et al. (2017); Riedinger et al. (2018). Each subsystem is a single-mode cavity with mode operators and that satisfy canonical commutation relations, . We now consider an arbitrary two-qubit bath state in the two-dimensional subspace spanned by ,

 |ψE⟩=bgg|gg⟩+bee|ee⟩. (19)

That is, the coefficients in Eq. (7) take values , while . For simplicity each cavity couples to its respective stream of bath qubits with the same rate () via an interaction that exchanges excitations, which corresponds to . Thus the ME in Eq. (9) has and

 ^L1 = √γ(bgg^a1+bee^a†2), (20a) ^L2 = √γ(bgg^a2+bee^a†1), (20b)

reminiscent of two-mode squeezing transformation.

The connection to two-mode squeezing can be made explicit under certain conditions. When we can define a strictly positive effective rate given by the population difference,

 Γ\coloneqqγ(|bgg|2−|bee|2). (21)

Then, the jump operators can be rewritten as

 ^L1 = √Γ[cosh(r)^a1+eiϑsinh(r)^a†2], (22a) ^L2 = √Γ[cosh(r)^a2+eiϑsinh(r)^a†1], (22b)

where the squeezing amplitude is related to the coefficients via the relations,

 cosh(r)=|bgg|√|bgg|2−|bee|2. (23)

The squeezing angle is given by the phase of relative to . Preparing the qubit bath in a maximally entangled state, , is the limit of infinite squeezing, . The jump operators are explicitly given by a two-mode squeezing transformation on the cavity annihilation operators, , where the unitary, two-mode squeezing operator is Wang et al. (2007); Schnabel (2017)

 ^S=eζ∗^a1^a2−ζ^a†1^a†2, (24)

with complex squeezing parameter . The transformation is in the Schrödinger-picture because the jump operators are nullifiers Gu et al. (2009); Menicucci et al. (2011) of the two-mode squeezed vacuum state. This can be seen directly by transforming the action of the annihilation operators,

 ^aℓ|0⟩2=0→^S^aℓ^S†^S|0⟩2∝^Lℓ|r,ϑ⟩2=0, (25)

where the two-mode squeezed vacuum state is

 |r,ϑ⟩2\coloneqq^S|0⟩2. (26)

Thus, the steady state of the ME is a pure, Gaussian, two-mode squeezed vacuum state with squeezing that depends on the qubit-bath coefficients. The dynamical preparation of from a two-mode vacuum state is shown for various values of the squeezing parameter in Fig. 2(a). The state is dissipatively cooled via interaction with the two-qubit bath towards the steady state, details can be found in Appendix B. We quantify the approach to with the Uhlmann-Jozsa fidelity, which can be calculated from the covariance matrix. As is increased, the effective rate decreases according to Eq. (21). That the time to approach the steady state increases exponentially with the squeezing parameter —see Fig. 2(b)—is unsurprising, since more highly squeezed states contain more energy.

In the opposite regime, where , there is no unitary Bogoliubov transformation that transforms the operators while maintaining the canonical commutation relations. Nevertheless, the jump operators in Eq. (20) may be written similarly to Eq. (22) with the roles of and reversed. Because , the jump operators contain a larger proportion of than , and the incoming two-qubit environment is more likely to transfer energy to the subsystems than to remove it. In this case the ME serves as an incoherent amplifier. In the following section we will investigate this parameter regime as well as the “exceptional points” where the bath is prepared in a maximally entangled state such as a Bell state.

Note that for the case with coefficients does not lead to the same dynamics even though the amount of entanglement in the bath can be the same. In this case and the resulting master equation just has dissipative terms not related to squeezing.

## V Two remote two-level atoms interacting with a Bell-state bath

In this section we investigate the repeated interaction between two remote two-level subsystems and a stream of maximally entangled bath qubits. The subsystems are taken to be identical, each described by a bare Hamiltonian , where is resonant with the bath qubit frequency, . The interaction between each subsystem and its bath qubit is an excitation exchange described by a lowering operator, . To avoid confusion, we henceforth refer to the bath as qubits and the subsystems as atoms, with the understanding that the bath qubits could indeed be physically manifested as a stream of entangled atoms.

### v.1 Bath qubits in a pure Bell state

We first consider a maximally entangled two-qubit bath state in the subspace ,

 |ψE⟩=1√2(|ee⟩+eiϕ|gg⟩), (27)

such that the Bell states and are given by , respectively. That is, the coefficients in Eq. (7) take values and while . For simplicity we set corresponding to decay rate , which yields the ME:

 ˙^ρ(t)= γ2(D[^σ1+eiϕ^σ†2]^ρ+D[eiϕ^σ†1+^σ2]^ρ). (28)

The joint state of the atomic subsystems is initialized in the arbitrary state

 ^ρ0=∑j,k∑j′,k′ρ0jk,j′k′|jk⟩⟨j′k′|, (29)

where the sums run over .

#### v.1.1 Bath qubits in |Φ+E⟩

First, we consider the case where the bath qubits are prepared in the Bell state, , given by in Eq. (27). In the long-time limit, , the steady state of the two-atom system is given by

 ^ρss=⎛⎜ ⎜ ⎜ ⎜⎝ρssee,ee00ρssee,gg0ρsseg,eg0000ρssge,ge0ρss∗ee,gg00ρssgg,gg⎞⎟ ⎟ ⎟ ⎟⎠, (30)

where the steady-state matrix elements are related to the initial-state matrix elements by

 ρssee,ee = 13[ρ0ee,ee+ρ0gg,gg−ρ0ee,gg+12ρ0eg,eg+12ρ0ge,ge], ρsseg,eg = 13[12ρ0ee,ee+12ρ0gg,gg+ρ0ee,gg+ρ0eg,eg+ρ0ge,ge], ρssge,ge = ρsseg,eg, ρssgg,gg = ρssee,ee, ρssee,gg = −16[ρ0ee,ee+ρ0gg,gg−4ρ0ee,gg−ρ0eg,eg−ρ0ge,ge].

Thus the atomic steady state is not unique; rather, the ME has an invariant subspace. A particular steady-state within this invariant subspace depends on the initial state Albert and Jiang (2014).

The steady state of the joint atomic system can exhibit non-classical correlations, identified by the negative partial transpose criterion Peres (1996). The partial transpose matrix , partitioned with respect to the subsystems, takes the following form,

 ^ρPTss=⎛⎜ ⎜ ⎜ ⎜⎝ρssee,ee0000ρsseg,egρssee,gg00ρssgg,eeρssge,ge0000ρssgg,gg⎞⎟ ⎟ ⎟ ⎟⎠. (31)

Negativity in the spectrum of guarantees the presence of entanglement. We quantify the entanglement by the logarithmic negativity Plenio (2005),

 LN(^ρ):=log2(Tr[√^ρ†PT^ρPT]), (32)

where with the minimum value corresponding to separable states and the maximum value to maximally entangled states.

A particular steady state of interest is when the atomic system is prepared in the Bell state , whence it does not undergo any decoherence via interaction with the qubit bath, since it is already at the steady state. However, when the initial atomic state is the steady state is the mixture,

 ^ρss=13(|Φ+⟩⟨Φ+|+|eg⟩⟨eg|+|ge⟩⟨ge|), (33)

with a positive partial transpose matrix, i.e. it has purely positive eigenvalues with boldface indicating degeneracy of order 3. A positive partial transpose is a sufficient condition for separability of a two-qubit system.

An interesting scenario is an initial joint atomic state, , that is a weighted sum of the two Bell states,

 |ψ0(θ)⟩=sinθ|Φ+⟩+cosθ|Φ−⟩, (34)

for . Beyond this range in , the entanglement behavior repeats. Above we found that the maximally entangled state, , is a steady state of the ME As deviates from zero, the contribution from the antisymmetric Bell state diminishes. Beyond the critical point () the effect of the symmetric Bell state is dominant and the entanglement vanishes, . At the critical point the steady state is a two-atom Werner state Werner (1989),

 ^ρss(θc)=13|Φ−⟩⟨Φ−|+16^1S, (35)

which is separable. In fact, even after this point the steady state is separable all the way to , at which point it is given by Eq. (33). The behavior is illustrated in Fig. 3(a) where we plot logarithmic negativity as a function of (dashed blue curve). Comparison with the initial logarithmic negativity (light gray curve) reveals that the Bell state environment preserves the initial system entanglement up until despite the fact that the state becomes mixed, Fig. 3(b).

#### v.1.2 Bath qubits in |Φ−E⟩, |Ψ+E⟩, or |Ψ−E⟩

The situation where the bath is prepared in proceeds similarly. In this case, the steady state has the same form as Eq. (30) with the following substitutions:

 ρ0ee,gg→−ρ0ee,gg,andρssee,gg→−ρssee,gg. (36)

When the initial atomic state is parameterized as in Eq. (34), the steady states for the extremal cases, , are just as in Eq. (33) and Eq. (35), respectively, with the roles of swapped. By varying in the interval , the antisymmetric Bell state is dominant, the result of which is that the systems remain separable. For the range the symmetric Bell state, is dominant in the initial atomic state, and the ME preserves the initial entanglement. For the initial state, , is a steady state.

Lastly, when the qubit environment is prepared in the other subspace ,

 |ψE⟩=1√2(|eg⟩+eiϕ|ge⟩), (37)

with Bell states and given by , the steady states have an analogous form. Table 1 summarizes the results for comparison.

### v.2 Bath qubits in a non-maximally entangled state

From the above analysis it may be inferred that the existence of coherences in the environment is merely a necessary condition for entangling the subsystems, even though the cross terms in the master equation, Eq. (28), might suggest otherwise, i.e. the generation of quantum correlation among systems. We found that Bell-state baths generate a ME without a unique steady state and, depending on the initial atomic state, atomic entanglement is either preserved or destroyed, but never created. Here we show that when the bath qubits are prepared in a non-maximally entangled state that can be arbitrarily close to a Bell state, the ME has a unique, entangled steady state. Consider the bath in the following state that slightly deviates from a maximally entangled state,

 |ψE(ϕ,ϵ)⟩=1√2+ϵ(|ee⟩+eiϕ