Creation of Two-Particle Entanglement in Open Macroscopic Quantum Systems

We consider an open quantum system of N not directly interacting spins (qubits) in contact with both local and collective thermal environments. The qubit-environment interactions are energy conserving. We trace out the variables of the thermal environments and N-2 qubits to obtain the time-dependent reduced density matrix for two arbitrary qubits. We numerically simulate the reduced dynamics and the creation of entanglement (concurrence) as a function of the parameters of the thermal environments and the number of qubits, N. Our results demonstrate that the two-qubit entanglement generally decreases as N increases. We show analytically that in the limit N tending to infinity, no entanglement can be created. This indicates that collective thermal environments cannot create two-qubit entanglement when many qubits are located within a region of the size of the environment coherence length. We discuss possible applications of our approach to the development of a new quantum characterization of noisy environments.


Introduction
In open many-body systems, such as solid-state and biological ones, quantum behavior reveals itself in many ways. Often the quantitative parameter used to measure "quantumness" (possibly even of macroscopic order) is entanglement. The presence of entanglement implies that the wave function (or the reduced density matrix) cannot be represented as a product of the corresponding objects for the individual qubits. It is important to note that to produce and to measure entanglement in such systems, one does not necessarily need to know much detail about the system, possibly not even its Hamiltonian [1]. The questions then are: How useful is entanglement as a measure of quatumness, what can it add to our knowledge of system properties and behavior, and how can it be utilized? Indeed, just knowing that the system is entangled (knowledge of complicated quantum behavior) is not sufficient to imply that its quantum properties are useful for specific applications. Fortunately, however, in some situations entanglement provides very useful properties, including additional exponential resources for quantum computation [2] and a possible enhancement of photosynthesis in bio-systems [3]. (See also references therein.) Entanglement could be produced by direct interaction between the qubits. This interaction should be of a "conditional" nature, mixing an initial product state (disentangled) in such a way that the final state becomes correlated in a quantum way. Many aspects of entanglement creation are widely discussed in the literature. (See, for example, [4][5][6][7][8][9] and references therein.) More recently, interest appeared in the possibility to create entanglement in the absence of direct interactions between qubits (or when the latter are very small). Entanglement can then still be created merely by the indirect interaction of non-interacting qubits through a collective thermal bath. In [4] this situation was considered for a model of two non-interacting spins 1/2 (qubits) interacting only with a collective thermal bosonic environment. It was demonstrated numerically in [4] that for some initially unentangled two-qubit states, and under some conditions on the thermal bath, measurable entanglement between the two qubits is created for intermediate times. The model of [4] is energy-conserving, ignoring relaxation processes for the qubits, and including only the effects of decoherence. In [9] these results were extended to a more general model having (i) both local and collective thermal environments (at the same temperature) and (ii) energy conserving and energy exchange interactions between qubits and their environments. The conditions for entanglement creation were discussed and analyzed numerically in [9]. It was concluded that, in spite of the competition between the local thermal environments (which destroy entanglement) and the collective thermal environment (tending to create entanglement), the creation of measurable entanglement can be realized for some finite times. In both papers, [4] and [9], only two qubits are analyzed, hence no direct connection to many-body systems was made. As was recently shown in [10], the presence of a large number of indirectly interacting qubits interacting only with their common collective thermal environment could significantly modify the effective single-qubit characteristics including their relaxation and decoherence rates.
In the present paper, we consider a model of N not directly interacting spins 1/2 (qubits) placed in a constant effective magnetic field (oriented in the z-direction). The qubits interact with both local and collective thermal environments (all at the same temperature). The collective interaction introduces an indirect qubit interaction. In the total density matrix of all qubits and environments, we trace over the variables of the environments and N − 2 qubits. This gives us the time-dependent reduced density matrix for two arbitrary qubits. In the z-representation it is represented by a timedependent 4 × 4 matrix. It is important to notice that the matrix elements, [ρ t ] n,m , n, m = 1, ..., 4, depend not only on the parameters of the thermal environments but also on the total number of qubits, N . We study numerically the concurrence C(t) of the reduced two-qubit density matrix, and its dependence on the parameters of the system and on N . To realize and study this situation in an experiment, one must have access to the two selected qubits (such as their particular frequencies), in order to manipulate them and prepare the initial state. Our main result is that the amplitude of concurrence, C max (t), generally decreases as N increases. This means that one should not expect that the collective thermal environment can create by itself measurable entanglement even of two-qubits, in the presence of many other qubits within a range of the collective environment coherence length.

Outline of main results
The initial state of the entire system is disentangled, a product state in which each of the N spins is in a state ρ j , j = 1, . . . , N , all local reservoir states are thermal and so is that of the collective reservoir (at a fixed common temperature 1 ).
Analytic results.
• Explicit dynamics. As the spins interact with the reservoirs via energy-conserving couplings only, the reduced two-spin dynamics can be calculated explicitly, see Proposition 2.1 below. Consequences of the energy conservation are that populations, i.e., the diagonal density matrix elements, are time-independent, and that the off-diagonal elements evolve independently. As an example, we discuss here the dynamics of the (1,2) matrix element, The other matrix elements have similar behavior. Each factor on the r.h.s. has an interpretation: -[ρ 0 ] 12 is the initial condition of the matrix element in question. None of the other initial matrix elements are involved (energy conserving coupling); -e iω 2 t is the uncoupled dynamics (no interaction with environments); -e iκ 2 c S(t) is a dephasing factor with a time-dependent phase S(t) ≤ 0 becoming linear for large t (for the considered infra-red behavior |k| 1/2 of the coupling constants in three dimensions, see (2.9)); it represents a "Lamb shift" contribution to the real part of the effective energy; this term is generated by the collective reservoir, but it is independent of the presence of the N − 2 traced-out spins (the term would be the same if only two spins were coupled to the reservoirs); -e −κ 2 c Γc(t) is a decaying factor with time-dependent decay rates, Γ(t) ≥ 0, becoming linear for large t (see (2.10)) both the local and collective reservoirs contribute; however, the term is independent of the N − 2 traced-out spins (again, it would be the same if only two spins were coupled to the reservoirs); -P N (t) is a product of N − 2 oscillating terms encoding the effect of all the tracedout spins (see (2.8)). It is important to notice that P N (t) only depends on the diagonal density matrix elements of the initial states of the N −2 traced-out spins. 2 Consequently, the two-qubit state does not depend on the initial off-diagonal density matrix elements of the N − 2 traced-out "background" spins. Typically, we expect those spins to be initially in (close to) equilibrium, corresponding to vanishing off-diagonals.
Some general properties of P N (t) can be explained easily for the case in which all N − 2 spins are initially in the high temperature equilibrium state 1 c S(t))] N −2 and its magnitude oscillates between zero and one. Due to the large power (N − 2), the peaks of the function |P N (t)|, centered at the discrete times t p satisfying S(t p ) ∈ π But the density matrix becomes very simple if P N (t) = 0, because many entries vanish (c.f. Proposition (2.1)) and the corresponding concurrence is zero. It follows that in the large N limit, concurrence is zero for all times (except possibly for some isolated instances, t p ).
• N -dependent scaling of the interaction κ c . The above analysis suggests that one cannot generate two-spin entanglement for N large at fixed interaction strength κ c . However, the width of the peaked function P N (t) which is of order . Hence we consider a N -dependent scaling of the coupling, replacing κ c by κ c /N η , for some η > 0. According to the above discussion, the borderline case is η = 1/4.
Starting from the explicit expressions (Proposition 2.1) and using the scaling κ c /N η , we calculate the limit N → ∞ of ρ t , for t ∈ R fixed. The analytic expressions we obtain for 0 < η < 1/4 and η > 1/4 show that the limiting dynamics does not create entanglement, for any time t. While we are able to obtain explicit expressions for concurrence in the regime of N → ∞, we are not so for N finite. 3 However, since no entanglement is generated in the limiting case, N → ∞, but we know entanglement is created for N = 2 (see e.g. [4,9]), we expect that entanglement creation decays with increasing N . We study this decay numerically.

Numerical results.
We introduce ν c , the highest frequency at which spin-reservoir interactions occur and call it the cutoff frequency. In the simulations, we take ν c of the order of the thermal frequency ν T = k B T /h. In the infra-red regime, our coupling is proportional to |k| (see after 2.10).
• For N = 2, concurrence creation is maximal if both spins start out in their high-temperature state 1 √ 2 (|+ + |− ), see Fig.1. Consequently, in the subsequent simulations, we take initial states of the two not traced-out qubits very close to this state, and we take the diagonals of the initial states of the N − 2 traced-out quibts to be constant 1/2 (remember that the off-diagonals of these qubits do not influence the dynamics at all). In Fig.4 we modify the initial state of the two not traced-out qubits and check that maximal concurrence is indeed obtained when both qubits are in the above state, even for large N .
• For general N , entanglement evolves according to a rescaled time t → κ 2 c ν c t, see Fig.2. This figure shows that a reduction of κ c diminishes the created concurrence in a moderate way. For instance decreasing κ c by a factor 10 only decreases concurrence by less than 1/3.
• In Fig.3 we show that the maximum of created concurrence decays with increasing N . For intermediate values of N (with the current parameters N ∼ 10 − 150) the decrease is exponential, for smaller and larger values of N , it is superexponential.
• In the same Fig.3c, we study the dependence of the maximal time, τ c , (before recurrence) at which the concurrence is not zero. We have found that this time decays exponentially in the number of spins, N , for sufficiently large N .
• Results on the rescaled model κ c → κ c /N η are shown in Fig.5. We find a decrease of maximal concurrence with increasing N for all η ≥ 0. The critical value, η = 1/4 (see analytic results above) divides the concurrence decay into two regimes. In the range, η > 1/4, the maximal concurrence decreases exponentially in N , for intermediate values of N (between 10 and 180), with a universal decay rate (i.e., not depending on η). For η < 1/4 the decay is superexponential and varies with η. We conclude that no scaling κ c → κ c /N η can compensate the decay of created concurrence for large N .

Model and reduced density matrix
The full Hamiltonian of the N noninteracting spins 1/2 coupled by energy conserving interactions to local and collective bosonic heat reservoirs is given by Below we use dimensionless variables and parameters. To do so, we introduce a characteristic frequency, ω 0 , typically of the order of spin transition frequency. The total Hamiltonian, energies of spin states, and temperature are measured in units ω 0 . The frequencies of spins, ω n > 0, bosonic excitations, ω(k) = c| k| (where c is the speed of light), the wave vectors of bosonic excitations are normalized by ω 0 /c, and all constants of interactions are measured in units of ω 0 . The dimensionless time is defined as ω 0 t.
In (2.1), (2.2), ω n > 0 is the frequency of spin n, and S z n denotes the S z of spin n. H R is the Hamiltonian of the bosonic collective reservoir, and H Rn is that same Hamiltonian for the n-th individual reservoir. For a squareintegrable form factor h(k), k ∈ R 3 , φ(h) is given by The real numbers, κ n and ν n , are coupling constants, measuring the strengths of the energy-conserving collective coupling and the energy-conserving local coupling, respectively.
Since the spins interact with the reservoirs only through energy-conserving channels, this model is exactly solvable. For simplicity of exposition, we take κ n = κ c for all n (collective) and ν n = ν ℓ for all n (local).
We also take, for simplicity, all local form factors equal (f ℓ ) and all collective ones also (f c ). Fix any pair of spins, and (re-)label their frequencies by ω 1 and ω 2 , see (2.1). We write the reduced density matrix, ρ t , of the two fixed spins as a 4 × 4 matrix [ρ t ] ij in the ordered energy basis The initial state of the spins is a product state of the form ρ S 1 ,0 ⊗ · · · ⊗ ρ S N ,0 , where with 0 ≤ p j ≤ 1 and |v j | 2 ≤ p j (1−p j ). The upper bound on the off-diagonal guarantees that the eigenvalues of ρ S j ,0 are non-negative. We introduce the quantities The integrals in (2.9) and (2.10) are made to converge introducing a suitable cutoff wavenumber, |k c |, or cut-off frequency, ν c = |k c |/2π. (Here we use dimensionless units.) For instance, for numerical simulations, we choose as form factor the function f c (k) = |k| χ |k|≤|kc| , where χ |k|≤|kc| = 1 if |k| ≤ |k c | and χ |k|≤|kc| = 0 otherwise. We also define P N (t) to be the same as P N (t), but with κ 2 c replaced by 2κ 2 c . With this notation, we have the following result.

Proposition 2.1 (Explicit dynamics of the reduced density matrix.) The evolution of the density matrix is given by
14) The proof of this proposition is a rather simple calculation. One can proceed as in [15] (proof of Proposition 7.4).
Remarks. 1. The effect of spins 3, . . . , N is contained entirely in the factors P N (t) and P N (t). They only depend on the initial populations p j , j = 3, . . . , N (see (2.8)), but not on the off-diagonals, v j . This is explained by the fact that when tracing over a single spin, j ≥ 3, we perform the operation Tr spinj U ρ S j ,0 V , where U, V are operators commuting with S z 2 (energy conserving interactions only!). Clearly the latter trace only involves the diagonal of ρ S j ,0 .
2. The oscillatory phases, e iωt , in (2.11)-(2.16) represent the free, uncoupled dynamics of the spins. Consider the modified two-spin density matrix ("interaction picture" dynamics of ρ t ). Because ρ ′ t and ρ t are related by conjugation of a unitary operator of the product form, e it(−ω 1 S z 1 ) ⊗ e it(−ω 2 S z 2 ) , the concurrences of ρ t and ρ ′ t are the same. In other words, when examining concurrence of ρ t , we may use formulas (2.11)-(2.16) with ω 1 = ω 2 = 0.

Variation of N -dependence
For homogeneous initial conditions, p j = p for j = 3, . . . , N , we have Unless p = 0, 1, |P N (t)| oscillates in t between its minimum value |2p − 1| N −2 (when cos(2κ 2 c S(t)) = −1) and its maximum value 1 (when cos(2κ 2 c S(t)) = 1). 5 The width of the oscillations becomes very narrow with increasing N . In the limit of large N , P N (t) is zero for all times, except for the discrete set of t ∈ R satisfying cos(2κ 2 c S(t)) = 1, in which case |P N (t)| = 1.
This implies that for large N , all off-diagonal density matrix elements of ρ t vanish with the exception of [ρ t ] 23 = e i(ω 1 −ω 2 )t e −2κ 2 ℓ Γ ℓ (t) (and [ρ t ] 32 , of course) for almost all values of t. This suppression of off-diagonals comes from the large number of particles and is mediated through the collective energy-conserving interaction. (For κ c = 0 we have P N (t) = 1.) In order to try to have a non-trivial dynamics for large N , one may scale the collective conserving coupling constant as Then (2.19) becomes (2.20) An expansion in large N yields Thus as N → ∞, (2.21) Remarks. 1. By replacing, in these limits, κ 2 c by 2κ 2 c , we obtain the corresponding limits for P N (t).

Asymptotic concurrence (N → ∞)
• 0 < η < 1/4 : The reduced two-spin density matrix (in the interaction picture) at time t is from which we obtain the concurrence Remember that |v j | 2 ≤ p j (1 − p j ). This shows that the N → ∞ asymptotic dynamics cannnot create entanglement at any time.
• η > 1/4 : Call the r.h.s. of (2.21) P ∞ (t) (a quantity still depending on N unless η = 1/2). By replacing κ 2 c by 2κ 2 c in (2.21) we obtain the limit of P N (t), which we call P ∞ (t). For η > 1/4 we have the relation P ∞ (t) = [P ∞ (t)] 2 . The reduced two-spin density matrix (in the interaction picture) at time t is, for N → ∞ The density matrix (2.24) is of the product form This shows that the N → ∞ asymptotic dynamics is factorizable and cannnot create entanglement at any time.

Numerical Results
• Let us first consider the case of two spins only, N = 2. In (2.11)-(2.16) we put for simplicity Γ ℓ = Γ c = Γ and regard Γ and S as two independent parameters. Taking both spins initially in the same state given by p, v, see (2.7), we examine the maximal concurrence, as a function of S and Γ, for arbitrary fixed values of p and v.
We find that for fixed p, v, the maximal concurrence is given at S = π/2, Γ = 0. Having such values fixed and plotting the concurrence as a function of p, v, the maximal concurrence is realized when p = v = 1/2, see Fig. 1, where a plot of the concurrence as a function of p, v is shown. Maximal generation of concurrence is thus obtained starting from pure state initial conditions 1 √ 2 (|+ + |− ) for each spin.
• Let us now consider the case of N spins. For concreteness we choose p = 1/2 for all spins (the traced-out ones and the two not traced-out ones). For the two not traced-out spins we take off-diagonals v = 0.48. (Then p is close to v which favors larger entanglement creation). Recall that the dynamics is independent of the off-diagonals of the N − 2 traced-out spins (i.e., we do not have to specify the v of the N − 2 traced-out spins).
As mentioned after (2.10), we choose the form factor, f c (k) = |k| χ |k|≤2πνc , with the cut-off frequency equal to the thermal frequency, ν c = ν T = k B T /h, at room temperature, T = 300 K.
In Fig. 2 we investigate the effects of an increase in the coupling parameter, κ c . The first effect is a time-shift for the concurrence evolution, described by a scaling, t → κ 2 c ν c t, see Fig. 2a, where ν c is the cut-off frequency. The second effect is a reduction of the maximal concurrence in a smooth way, see Fig. 2b. As one can see, the effective decrease in amplitude for not too strong coupling strength, κ c , is relatively small. For instance, changing κ c for N = 2 by one order of magnitude from 0.04 to 0.4, changes the amplitude by only 27%. The percentage change is almost the same for larger N values, see Fig. 2b.
In Fig. 3 we show that the creation of concurrence decreases with the number of spins. In graph a) we plot the concurrence as a function of (rescaled) time for various values of N = 2, . . . , 32. As one can see, the same time rescaling is also valid for N > 2. Moreover, the maximum concurrence created, C max , reported in graph b), decreases exponentially in N in the range 10 < N < 150 and faster than exponentially outside this range. For larger N , the concurrence decays superexponentially in N . For N exceeding 200, the concurrence becomes too small to be significant (of the order 10 −4 ).
It is also interesting to note that the graph of concurrence shows collapses and revivals and that the revival times for N > 2 are always less than the revival time for N = 2. It is also interesting to consider how the collapse time, τ c , defined as the first time at which concurrence drops abruptly to zero, depends on the number of spins, N . This study has been reported in Fig. 3c  • One can also vary the initial conditions for the spins by choosing independent p 1,2 and v 1,2 , while all other spins have the same value p j = 1/2, j = 3, . . . , N (their off-diagonals v j do not influence the dynamics at all).
In order to simplify the problem, we also set p 1 = v 1 and p 2 = v 2 and consider the maximal concurrence as a function of two independent parameters, p 1 and p 2 , only.
An example of the 3D plot obtained is reported in Fig. 4 : the maximal concurrence is realized for p 1 = p 2 = v 1 = v 2 = 1/2, independently of N (in the picture N = 40, but similar plots are obtained for other values of N ).
• The numerical analysis of the rescaled model with κ c replaced by κ c /N η shows that the concurrence is always a decreasing function of N and that the maximum of the created concurrence is a universal function of the number of spins N , independent of η for η > 1/4.
Results are shown in Fig. 5, where the dashed line is the best exponential fit exp(−aN ), with a = 0.0177 ± 0.0003, the best fitting value, for the cases η > 1/4. The same figure shows that when η ≤ 1/4 the decay is superexponential and no universality occurs.
This suggests that no power law scaling with N of the coupling strength can compensate the rapid decay of concurrence with the number of spins. Qualitatively similar results, not reported here, can be obtained by changing the ratio between the thermal and cut-frequency in the range (0.5, 4).

Conclusion
We have analyzed the two-qubit (two effective spin) entanglement in an N -qubit open system, where individual qubits interacting through a collective thermal environment.
We have demonstrated that concurrence (a measure of two-qubit entanglement) quickly decays with increasing number of surrounding qubits. It follows from our consideration that for creation of entanglement by implementing our approach, one has to use a small number of qubits collectively interacting with the thermal environment, preferably only two qubits. In this paper, we consider a quantum noise generated by a thermal bosonic environment, interacting in a purely dephasing way with qubits, but no Markovian approximation has been made. In contrast, it has been shown [11][12][13] that an N -qubit register can be driven to entangled target states (within a wide class of states), under an appropriately engineered Markovian dissipative dynamics in Lindblad form, and that the convergence rate to the target state is independent of N . This entirely different noise has thus a completely different effect on entanglement.
The results of this paper are important for a better understanding and a characterization of the collective thermal environment and its ability to create and destroy entanglement in many-particle open quantum systems. Any quantum information processor could be used for the study of entanglement induced by a thermal environment, as for example two (or more) ions in an ion trap quantum computer [2] or two (or more) superconducting qubits [14], or an NMR quantum information processor [2] operating at room temperature. In the latter case, the experiment is performed with an ensemble of many identical molecules, each having a similar number of weakly interacting nuclear spins. One can study in this system, for example, how the collective vibrational modes of a molecule (thermal environment) influence the two-spin "effective entanglement" in the ensemble of molecules, depending on the number of spins in a molecule. Important and relevant systems are represented by photosynthetic complexes in which quantum effects for exciton dynamics and primary charge separation have been recently discovered [16][17][18][19]. Here, the effective spins are exciton, or electron sites (pigments) and the thermal environment represents the vibrational protein modes. Here, N = 2 spins is considered.