# Creation of Two-Particle Entanglement in Open Macroscopic Quantum Systems

Authors: M. Merkli; G. P. Berman; F. Borgonovi; V. I. Tsifrinovich
Journal: Advances in Mathematical Physics
 (2012)

---

## Abstract

We consider an open quantum system ofN 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→∞, 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 relevance of our consideration
to recent quantum information devices and biosystems.


---

## Body

## 1. 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, e.g., [4–11] 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 noninteracting qubits through a collective thermal bath. In [4] this situation was considered for a model of two noninteracting 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 [11] 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 [11]. 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, 11], only two qubits are analyzed hence no direct connection to many-body systems was made. As was recently shown in [12], 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 ofN 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 time-dependent 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, Cmax(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 the range of the collective environment coherence length.


### 1.1. Outline of Main Results

The initial state of the entire system is disentangled, a product state in which each of theN 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 (A generalization to a nonequilibrium situation where each local and the collective reservoir have different individual temperatures is immediate.).


#### 1.1.1. Analytic Results

(i) Explicit Dynamics
As the spins interact with the reservoirs via energy-conserving couplings only, the reduced two-spin dynamics can be calculated explicitly; see Proposition2.1. Consequences of the energy conservation are that populations, that is, 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:
(1.1)[ρt]12=[ρ0]12eiω2teiϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t).
The other matrix elements have similar behavior. Each factor on the r.h.s. has an interpretation.(i)
[ρ0]12 is the initial condition of the matrix element in question. None of the other initial matrix elements are involved (energy-conserving coupling).(ii)
eiω2t is the uncoupled dynamics (no interaction with environments);(iii)
eiϰc2S(t) is a dephasing factor with a time-dependent phase S(t)≤0 becoming linear for large t (for the considered infrared behavior |k|1/2 of the coupling constants in three dimensions, see (2.10)); 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).(iv)
e-ϰℓ2Γℓ(t)-ϰc2Γc(t) is a decaying factor with time-dependent decay rates, Γ(t)≥0, becoming linear for large t (see (2.11)). 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);(v)
PN(t) is a product of N-2 oscillating terms encoding the effect of all the traced-out spins (see (2.9)). It is important to notice that PN(t) only depends on the diagonal density matrix elements of the initial states of theN-2traced-out spins. (This is so since the dynamics is energy-conserving, and tracing out any of the spins involves only the diagonal of the initial (time zero) density matrix. (See also Remark 2.2(1) after Proposition 2.1.)). 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 ofPN(t) can be explained easily for the case in which all N-2 spins are initially in the high-temperature equilibrium state (1/2)(|+〉+|-〉). Then PN(t)=[cos(ϰc2S(t))]N-2 and its magnitude oscillate between zero and one. Due to the large power (N-2), the peaks of the function |PN(t)|, centered at the discrete times tp satisfying S(tp)∈(π/ϰc2)ℤ, are of very narrow width O(1/(ϰc2N)) for large N. Consequently, in the limit N→∞, with ϰc held fixed, |PN(t)| is zero for all t except for t=tp, where |PN(tp)|=1. But the density matrix becomes very simple if PN(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, tp).(ii)N-Dependent Scaling of the Interaction ϰc
The above analysis suggests that one cannot generate two-spin entanglement forN large at fixed interaction strength ϰc. However, the width of the peaked function PN(t) which is of order O(1/(ϰc2N)) becomes appreciable if ϰc≳N1/4. Hence we consider an 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 (Proposition2.1) and using the scaling ϰc/Nη, we calculate the limit N→∞ of ρt, for t∈ℝ 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. (The reduced density matrix is given explicitly for all N and all t, but calculating from it explicitly the concurrence is more difficult.) However, since no entanglement is generated in the limiting case, N→∞, but we know entanglement is created for N=2 (see, e.g., [4, 11]), we expect that entanglement creation decays with increasing N. We study this decay numerically.


#### 1.1.2.   Numerical Results

We introduceνc, the highest frequency at which spin-reservoir interactions occur and call it the cut-off frequency. In the simulations, we take νc of the order of the thermal frequency νT=kBT/h. In the infra-red regime, our coupling is proportional to |k| (see after (2.11)).(
i
) For N=2, concurrence creation is maximal if both spins start out in their high-temperature state (1/2)(|+〉+|-〉); see Figure 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 Figure 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.Figure 1
Maximal concurrence as a function ofp and v, for fixed Γ(t)=0 ans S(t)=π/2. Here, N=2 spins is considered.(
ii
) For general N, entanglement evolves according to a rescaled time t↦ϰc2νct, see Figure 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.(a) Concurrence as a function of the rescaled time, for differentϰc and fixed ϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor fc(k)=|k|. Other data are p=0.5,v=0.48, N=4. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the coupling strength, ϰc, for different N values as indicated in the legend. Other values are the same as in (a).
(a)
(b)(
iii
) In Figure 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.(a) Concurrence as a function of the rescaled time, for fixedϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor, fc(k)=|k|. Other data are p=0.5, v=0.48, ϰc=0.05. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the number of spins N. (c) Plot of the collapse time, τc, as a function of the number of spins, N. The dashed line is the best fit exp(-αN) with α=0.0838±0.0002.
(a)
(b)
(c)Figure 4
Maximal concurrence as a function of the independent parameters,p1=v1 and p2=v2. As one can see the maximal concurrence is realized at the external corner, that is  p1=p2=1/2. Here is ϵ=1,  N=40 and p=1/2 for all other spins.(
iv
) In the same Figure 3(c), 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.(
v
) Results on the rescaled model ϰc↦ϰc/Nη are shown in Figure 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.Figure 5
Maximal concurrence as a function of the number of spins, for different power law scaling, as indicated in the legend. Here isfc(k)=|k|, p=0.5, v=0.48, ϵ=νc/νT=1, T=300 K, ϰc=0.2. The dashed line indicates a fitting exponential for the cases η>1/4. The solid curve indicates the case, η=1/4.


## 1.1. Outline of Main Results

The initial state of the entire system is disentangled, a product state in which each of theN 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 (A generalization to a nonequilibrium situation where each local and the collective reservoir have different individual temperatures is immediate.).


### 1.1.1. Analytic Results

(i) Explicit Dynamics
As the spins interact with the reservoirs via energy-conserving couplings only, the reduced two-spin dynamics can be calculated explicitly; see Proposition2.1. Consequences of the energy conservation are that populations, that is, 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:
(1.1)[ρt]12=[ρ0]12eiω2teiϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t).
The other matrix elements have similar behavior. Each factor on the r.h.s. has an interpretation.(i)
[ρ0]12 is the initial condition of the matrix element in question. None of the other initial matrix elements are involved (energy-conserving coupling).(ii)
eiω2t is the uncoupled dynamics (no interaction with environments);(iii)
eiϰc2S(t) is a dephasing factor with a time-dependent phase S(t)≤0 becoming linear for large t (for the considered infrared behavior |k|1/2 of the coupling constants in three dimensions, see (2.10)); 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).(iv)
e-ϰℓ2Γℓ(t)-ϰc2Γc(t) is a decaying factor with time-dependent decay rates, Γ(t)≥0, becoming linear for large t (see (2.11)). 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);(v)
PN(t) is a product of N-2 oscillating terms encoding the effect of all the traced-out spins (see (2.9)). It is important to notice that PN(t) only depends on the diagonal density matrix elements of the initial states of theN-2traced-out spins. (This is so since the dynamics is energy-conserving, and tracing out any of the spins involves only the diagonal of the initial (time zero) density matrix. (See also Remark 2.2(1) after Proposition 2.1.)). 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 ofPN(t) can be explained easily for the case in which all N-2 spins are initially in the high-temperature equilibrium state (1/2)(|+〉+|-〉). Then PN(t)=[cos(ϰc2S(t))]N-2 and its magnitude oscillate between zero and one. Due to the large power (N-2), the peaks of the function |PN(t)|, centered at the discrete times tp satisfying S(tp)∈(π/ϰc2)ℤ, are of very narrow width O(1/(ϰc2N)) for large N. Consequently, in the limit N→∞, with ϰc held fixed, |PN(t)| is zero for all t except for t=tp, where |PN(tp)|=1. But the density matrix becomes very simple if PN(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, tp).(ii)N-Dependent Scaling of the Interaction ϰc
The above analysis suggests that one cannot generate two-spin entanglement forN large at fixed interaction strength ϰc. However, the width of the peaked function PN(t) which is of order O(1/(ϰc2N)) becomes appreciable if ϰc≳N1/4. Hence we consider an 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 (Proposition2.1) and using the scaling ϰc/Nη, we calculate the limit N→∞ of ρt, for t∈ℝ 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. (The reduced density matrix is given explicitly for all N and all t, but calculating from it explicitly the concurrence is more difficult.) However, since no entanglement is generated in the limiting case, N→∞, but we know entanglement is created for N=2 (see, e.g., [4, 11]), we expect that entanglement creation decays with increasing N. We study this decay numerically.


### 1.1.2.   Numerical Results

We introduceνc, the highest frequency at which spin-reservoir interactions occur and call it the cut-off frequency. In the simulations, we take νc of the order of the thermal frequency νT=kBT/h. In the infra-red regime, our coupling is proportional to |k| (see after (2.11)).(
i
) For N=2, concurrence creation is maximal if both spins start out in their high-temperature state (1/2)(|+〉+|-〉); see Figure 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 Figure 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.Figure 1
Maximal concurrence as a function ofp and v, for fixed Γ(t)=0 ans S(t)=π/2. Here, N=2 spins is considered.(
ii
) For general N, entanglement evolves according to a rescaled time t↦ϰc2νct, see Figure 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.(a) Concurrence as a function of the rescaled time, for differentϰc and fixed ϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor fc(k)=|k|. Other data are p=0.5,v=0.48, N=4. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the coupling strength, ϰc, for different N values as indicated in the legend. Other values are the same as in (a).
(a)
(b)(
iii
) In Figure 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.(a) Concurrence as a function of the rescaled time, for fixedϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor, fc(k)=|k|. Other data are p=0.5, v=0.48, ϰc=0.05. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the number of spins N. (c) Plot of the collapse time, τc, as a function of the number of spins, N. The dashed line is the best fit exp(-αN) with α=0.0838±0.0002.
(a)
(b)
(c)Figure 4
Maximal concurrence as a function of the independent parameters,p1=v1 and p2=v2. As one can see the maximal concurrence is realized at the external corner, that is  p1=p2=1/2. Here is ϵ=1,  N=40 and p=1/2 for all other spins.(
iv
) In the same Figure 3(c), 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.(
v
) Results on the rescaled model ϰc↦ϰc/Nη are shown in Figure 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.Figure 5
Maximal concurrence as a function of the number of spins, for different power law scaling, as indicated in the legend. Here isfc(k)=|k|, p=0.5, v=0.48, ϵ=νc/νT=1, T=300 K, ϰc=0.2. The dashed line indicates a fitting exponential for the cases η>1/4. The solid curve indicates the case, η=1/4.


## 1.1.1. Analytic Results

(i) Explicit Dynamics
As the spins interact with the reservoirs via energy-conserving couplings only, the reduced two-spin dynamics can be calculated explicitly; see Proposition2.1. Consequences of the energy conservation are that populations, that is, 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:
(1.1)[ρt]12=[ρ0]12eiω2teiϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t).
The other matrix elements have similar behavior. Each factor on the r.h.s. has an interpretation.(i)
[ρ0]12 is the initial condition of the matrix element in question. None of the other initial matrix elements are involved (energy-conserving coupling).(ii)
eiω2t is the uncoupled dynamics (no interaction with environments);(iii)
eiϰc2S(t) is a dephasing factor with a time-dependent phase S(t)≤0 becoming linear for large t (for the considered infrared behavior |k|1/2 of the coupling constants in three dimensions, see (2.10)); 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).(iv)
e-ϰℓ2Γℓ(t)-ϰc2Γc(t) is a decaying factor with time-dependent decay rates, Γ(t)≥0, becoming linear for large t (see (2.11)). 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);(v)
PN(t) is a product of N-2 oscillating terms encoding the effect of all the traced-out spins (see (2.9)). It is important to notice that PN(t) only depends on the diagonal density matrix elements of the initial states of theN-2traced-out spins. (This is so since the dynamics is energy-conserving, and tracing out any of the spins involves only the diagonal of the initial (time zero) density matrix. (See also Remark 2.2(1) after Proposition 2.1.)). 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 ofPN(t) can be explained easily for the case in which all N-2 spins are initially in the high-temperature equilibrium state (1/2)(|+〉+|-〉). Then PN(t)=[cos(ϰc2S(t))]N-2 and its magnitude oscillate between zero and one. Due to the large power (N-2), the peaks of the function |PN(t)|, centered at the discrete times tp satisfying S(tp)∈(π/ϰc2)ℤ, are of very narrow width O(1/(ϰc2N)) for large N. Consequently, in the limit N→∞, with ϰc held fixed, |PN(t)| is zero for all t except for t=tp, where |PN(tp)|=1. But the density matrix becomes very simple if PN(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, tp).(ii)N-Dependent Scaling of the Interaction ϰc
The above analysis suggests that one cannot generate two-spin entanglement forN large at fixed interaction strength ϰc. However, the width of the peaked function PN(t) which is of order O(1/(ϰc2N)) becomes appreciable if ϰc≳N1/4. Hence we consider an 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 (Proposition2.1) and using the scaling ϰc/Nη, we calculate the limit N→∞ of ρt, for t∈ℝ 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. (The reduced density matrix is given explicitly for all N and all t, but calculating from it explicitly the concurrence is more difficult.) However, since no entanglement is generated in the limiting case, N→∞, but we know entanglement is created for N=2 (see, e.g., [4, 11]), we expect that entanglement creation decays with increasing N. We study this decay numerically.


## 1.1.2.   Numerical Results

We introduceνc, the highest frequency at which spin-reservoir interactions occur and call it the cut-off frequency. In the simulations, we take νc of the order of the thermal frequency νT=kBT/h. In the infra-red regime, our coupling is proportional to |k| (see after (2.11)).(
i
) For N=2, concurrence creation is maximal if both spins start out in their high-temperature state (1/2)(|+〉+|-〉); see Figure 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 Figure 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.Figure 1
Maximal concurrence as a function ofp and v, for fixed Γ(t)=0 ans S(t)=π/2. Here, N=2 spins is considered.(
ii
) For general N, entanglement evolves according to a rescaled time t↦ϰc2νct, see Figure 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.(a) Concurrence as a function of the rescaled time, for differentϰc and fixed ϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor fc(k)=|k|. Other data are p=0.5,v=0.48, N=4. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the coupling strength, ϰc, for different N values as indicated in the legend. Other values are the same as in (a).
(a)
(b)(
iii
) In Figure 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.(a) Concurrence as a function of the rescaled time, for fixedϵ=νc/νT=1, being respectively the cut-off and thermal frequency (νT=kBT/h with T=300K) and the form factor, fc(k)=|k|. Other data are p=0.5, v=0.48, ϰc=0.05. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the number of spins N. (c) Plot of the collapse time, τc, as a function of the number of spins, N. The dashed line is the best fit exp(-αN) with α=0.0838±0.0002.
(a)
(b)
(c)Figure 4
Maximal concurrence as a function of the independent parameters,p1=v1 and p2=v2. As one can see the maximal concurrence is realized at the external corner, that is  p1=p2=1/2. Here is ϵ=1,  N=40 and p=1/2 for all other spins.(
iv
) In the same Figure 3(c), 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.(
v
) Results on the rescaled model ϰc↦ϰc/Nη are shown in Figure 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.Figure 5
Maximal concurrence as a function of the number of spins, for different power law scaling, as indicated in the legend. Here isfc(k)=|k|, p=0.5, v=0.48, ϵ=νc/νT=1, T=300 K, ϰc=0.2. The dashed line indicates a fitting exponential for the cases η>1/4. The solid curve indicates the case, η=1/4.


## 2. Model and Reduced Density Matrix

The full Hamiltonian of theN noninteracting spins 1/2 coupled by energy conserving interactions to local and collective bosonic heat reservoirs is given by
(2.1)H=-ħ∑n=1N‍ωnSnz+∑n=1N‍HRn+HR(2.2)+∑n=1N‍ϰnSnz⊗ϕc(fc)+∑n=1N‍νnSnz⊗ϕn(fn).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 ω0t.In (2.1) and  (2.2), ωn>0 is the frequency of spin n,
(2.3)Sz=12[100-1],
and Snz denotes the Sz of spin n.  HR is the Hamiltonian of the bosonic collective reservoir,
(2.4)HR=∫R3‍|k|a*(k)a(k)d3k,
and HRn is that same Hamiltonian for the nth individual reservoir. For a square-integrable form factorh(k), k∈ℝ3, ϕ(h) is given by
(2.5)ϕ(h)=12∫R3‍{h(k)a*(k)+h(k)*a(k)}d3k.
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 isexactly solvable. For simplicity of exposition, we take
(2.6)ϰn=ϰc∀n(collective),νn=νl∀n(local).
We also take, for simplicity, all local form factors equal (fℓ) and all collective ones also (fc).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 (2.7)Φ1=φ1⊗φ1,Φ2=φ1⊗φ-1,Φ3=φ-1⊗φ1,Φ4=φ-1⊗φ-1,(Another equivalent notation is: Φ1=|++〉,Φ2=|+-〉,Φ3=|-+〉,Φ4=|--〉.) where Szφ±1=±(1/2)φ±1. For instance, [ρt]2,4=〈Φ2,ρtΦ4〉.Theinitial state of the spins is a product state of the form ρS1,0⊗⋯⊗ρSN,0, where
(2.8)ρSj,0=[pjvjvj*1-pj],
with 0≤pj≤1 and |vj|2≤pj(1-pj). The upper bound on the off-diagonal guarantees that the eigenvalues of ρSj,0 are nonnegative.We introduce the quantities:(2.9)PN(t)=∏j=3N‍[pjeiϰc2S(t)+(1-pj)e-iϰc2S(t)],(2.10)S(t)=-12∫R3‍|fc(k)|2|k|t-sin⁡(|k|t)|k|2d3k,(2.11)Γl,c(t)=∫R3‍|fl,c(k)|2coth⁡(β|k|2)sin⁡2(|k|t/2)|k|2d3k.The integrals in (2.10) and (2.11) are made to converge introducing a suitable cut-off wavenumber, |kc|, or cut-off frequency, νc=|kc|/2π. (Here we use dimensionless units.) For instance, for numerical simulations, we choose as form factor the function fc(k)=|k|χ|k|≤|kc|, where χ|k|≤|kc|=1 if |k|≤|kc| and   χ|k|≤|kc|=0 otherwise.We also defineP~N(t) to be the same as PN(t), but with ϰc2 replaced by 2ϰc2. 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(2.12)[ρt]12=[ρ0]12eiω2teiϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t),[ρt]13=[ρ0]13eiω1teiϰc2S(t)e-ϰc2Γc(t)PN(t),[ρt]14=[ρ0]14ei(ω1+ω2)te-2ϰl2Γl(t)-4ϰc2Γc(t)P~N(t),[ρt]23=[ρ0]23ei(ω1-ω2)te-2ϰl2Γl(t),[ρt]24=[ρ0]24eiω1te-iϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t),[ρt]34=[ρ0]34eiω2te-iϰc2S(t)e-ϰl2Γl(t)-ϰc2Γc(t)PN(t),
and the populations are constant, [ρt]jj=[ρ0]jj, for j=1,…,4 and t∈ℝ.The proof of this proposition is a rather simple calculation. One can proceed as in [13–15]  (proof of Proposition  7.4 in [13]).Remark 2.2.
(1) The effect of spins3,…,N is contained entirely in the factors PN(t) and P~N(t). They only depend on the initial populations pj, j=3,…,N (see (2.9)), but not on the off-diagonals, vj. This is explained by the fact that when tracing over a single spin, j≥3, we perform the operation TrspinjUρSj,0V, where U,V are operators commuting with S2z (energy-conserving interactions only!). Clearly the latter trace only involves the diagonal of ρSj,0.
(2) The oscillatory phases,eiωt, in (2.12) represent the free, uncoupled dynamics of the spins. Consider the modified two-spin density matrix:
(2.13)ρt′=eit(-ω1S1z-ω2S2z)ρte-it(-ω1S1z-ω2S2z),
(“interaction picture” dynamics of ρt). Because ρt′ and ρt are related by conjugation of a unitary operator of the product form, eit(-ω1S1z)⊗eit(-ω2S2z), the concurrences of ρt and ρt′ are the same. In other words, when examining concurrence of ρt, we may use formulas (2.12) with ω1=ω2=0.


### 2.1. Concurrence

Recall that the concurrence of the reduced density matrix,ρt, is unchanged when we pass to the interaction picture ρt′ (see the remark explaining (2.13)). In the basis (2.7), the evolution of ρt′ (2.13)  is given by (2.12) with ω1=ω2=0 and where the initial condition is (ρ0′=ρ0):
(2.14)[ρ0]=[p1p2p1v2v1p2v1v2p1v2*p1(1-p2)v1v2*v1(1-p2)v1*p2v1*v2(1-p1)p2(1-p1)v2v1*v2*v1*(1-p2)(1-p1)v2*(1-p1)(1-p2)].


#### 2.1.1. Variation ofN Dependence

For homogeneous initial conditions,pj=p for j=3,…,N, we have
(2.15)PN(t)=[peiϰc2S(t)+(1-p)e-iϰc2S(t)]N-2.
Unless p=0,1, |PN(t)| oscillates in t between its minimum value |2p-1|N-2 (when cos(2ϰc2S(t))=-1) and its maximum value 1 (when cos(2ϰc2S(t))=1) (We have |peiϰc2S(t)+(1-p)e-iϰc2S(t)|2=p2+2p(1-p)cos(2ϰc2S(t))+(1-p)2.) The width of the oscillations becomes very narrow with increasing N. In the limit of large N, PN(t) is zero for all times, except for the discrete set of t∈ℝ satisfying cos(2ϰc2S(t))=1, in which case |PN(t)|=1.This implies that, for largeN, all off-diagonal density matrix elements of ρt vanish with the exception of [ρt]23=ei(ω1-ω2)te-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 PN(t)=1.)In order to try to have a nontrivial dynamics for largeN, one may scale the collective conserving coupling constant as
(2.16)ϰc⟶ϰcNη,someη>0.
Then (2.15) becomes
(2.17)PN(t)=[peiϰc2S(t)/N2η+(1-p)e-iϰc2S(t)/N2η]N-2.
An expansion in large N yields
(2.18)PN(t)=e-iϰc2S(t)[1-2p]N1-2ηe-2ϰc4S2(t)N1-4η[p(1-p)+O(N-2η)].
Thus as N→∞,
(2.19)PN(t)⟶{e-iϰc2S(t)[1-2p]N1/2e-2ϰc4S2(t)p(1-p)N1-4η,0<η≤14(p≠0,1),e-iϰc2S(t)[1-2p]N1-2η,14<η.Remark 2.3.
(1) By replacing, in these limits,ϰc2 by 2ϰc2, we obtain the corresponding limits for P~N(t).
(2) For0<η<1/4,  PN(t) vanishes as N→∞.
(3) The rapid oscillating phases disappear ifp=1/2 (any η>0) or η=1/2 (any p).


#### 2.1.2. Asymptotic Concurrence (N→∞)

(
i
)
0
<
η
<
1
/
4
The reduced two-spin density matrix (in the interaction picture) at timet is
(2.20)[ρt′]=[p1p20000p1(1-p2)v1v2*e-2ϰl2Γl(t)00v1*v2e-2ϰl2Γl(t)(1-p1)p20000(1-p1)(1-p2)],
from which we obtain the concurrence
(2.21)C(ρt)=max⁡{0,-2[p1(1-p1)p2(1-p2)-|v1||v2|e-2ϰl2Γl(t)]}=0.
Remember that |vj|2≤pj(1-pj). This shows that the N→∞ asymptotic dynamics cannot create entanglement at any time.(
ii
)
η
>
1
/
4
Call the r.h.s. of (2.19) P∞(t) (a quantity still depending on N unless η=1/2). By replacing ϰc2 by 2ϰc2 in (2.19), 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→∞,
(2.22)[ρt′]=[p1p2p1v2DlP∞(t)v1p2DlP∞(t)v1v2Dl2P∞(t)2p1v2*DlP∞(t)*p1(1-p2)v1v2*Dl2v1(1-p2)DlP∞(t)v1*p2DlP∞(t)*v1*v2Dl2(1-p1)p2(1-p1)v2DlP∞(t)v1*v2*Dl2[P∞(t)2]*v1*(1-p2)DlP∞(t)*(1-p1)v2*DlP∞(t)*(1-p1)(1-p2)],
where
(2.23)Dl=Dl(t)=e-ϰl2Γl(t),P∞(t)=e-iϰc2S(t)[1-2p]N1-2η.
The density matrix (2.22) is of the product form:
(2.24)[p1v1Dl(t)P∞(t)v1*Dl(t)P∞(t)*1-p1]⊗[p2v2Dl(t)P∞(t)v2*Dl(t)P∞(t)*1-p2].
This shows that the N→∞ asymptotic dynamics is factorizable and cannot create entanglement at any time.


## 2.1. Concurrence

Recall that the concurrence of the reduced density matrix,ρt, is unchanged when we pass to the interaction picture ρt′ (see the remark explaining (2.13)). In the basis (2.7), the evolution of ρt′ (2.13)  is given by (2.12) with ω1=ω2=0 and where the initial condition is (ρ0′=ρ0):
(2.14)[ρ0]=[p1p2p1v2v1p2v1v2p1v2*p1(1-p2)v1v2*v1(1-p2)v1*p2v1*v2(1-p1)p2(1-p1)v2v1*v2*v1*(1-p2)(1-p1)v2*(1-p1)(1-p2)].


### 2.1.1. Variation ofN Dependence

For homogeneous initial conditions,pj=p for j=3,…,N, we have
(2.15)PN(t)=[peiϰc2S(t)+(1-p)e-iϰc2S(t)]N-2.
Unless p=0,1, |PN(t)| oscillates in t between its minimum value |2p-1|N-2 (when cos(2ϰc2S(t))=-1) and its maximum value 1 (when cos(2ϰc2S(t))=1) (We have |peiϰc2S(t)+(1-p)e-iϰc2S(t)|2=p2+2p(1-p)cos(2ϰc2S(t))+(1-p)2.) The width of the oscillations becomes very narrow with increasing N. In the limit of large N, PN(t) is zero for all times, except for the discrete set of t∈ℝ satisfying cos(2ϰc2S(t))=1, in which case |PN(t)|=1.This implies that, for largeN, all off-diagonal density matrix elements of ρt vanish with the exception of [ρt]23=ei(ω1-ω2)te-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 PN(t)=1.)In order to try to have a nontrivial dynamics for largeN, one may scale the collective conserving coupling constant as
(2.16)ϰc⟶ϰcNη,someη>0.
Then (2.15) becomes
(2.17)PN(t)=[peiϰc2S(t)/N2η+(1-p)e-iϰc2S(t)/N2η]N-2.
An expansion in large N yields
(2.18)PN(t)=e-iϰc2S(t)[1-2p]N1-2ηe-2ϰc4S2(t)N1-4η[p(1-p)+O(N-2η)].
Thus as N→∞,
(2.19)PN(t)⟶{e-iϰc2S(t)[1-2p]N1/2e-2ϰc4S2(t)p(1-p)N1-4η,0<η≤14(p≠0,1),e-iϰc2S(t)[1-2p]N1-2η,14<η.Remark 2.3.
(1) By replacing, in these limits,ϰc2 by 2ϰc2, we obtain the corresponding limits for P~N(t).
(2) For0<η<1/4,  PN(t) vanishes as N→∞.
(3) The rapid oscillating phases disappear ifp=1/2 (any η>0) or η=1/2 (any p).


### 2.1.2. Asymptotic Concurrence (N→∞)

(
i
)
0
<
η
<
1
/
4
The reduced two-spin density matrix (in the interaction picture) at timet is
(2.20)[ρt′]=[p1p20000p1(1-p2)v1v2*e-2ϰl2Γl(t)00v1*v2e-2ϰl2Γl(t)(1-p1)p20000(1-p1)(1-p2)],
from which we obtain the concurrence
(2.21)C(ρt)=max⁡{0,-2[p1(1-p1)p2(1-p2)-|v1||v2|e-2ϰl2Γl(t)]}=0.
Remember that |vj|2≤pj(1-pj). This shows that the N→∞ asymptotic dynamics cannot create entanglement at any time.(
ii
)
η
>
1
/
4
Call the r.h.s. of (2.19) P∞(t) (a quantity still depending on N unless η=1/2). By replacing ϰc2 by 2ϰc2 in (2.19), 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→∞,
(2.22)[ρt′]=[p1p2p1v2DlP∞(t)v1p2DlP∞(t)v1v2Dl2P∞(t)2p1v2*DlP∞(t)*p1(1-p2)v1v2*Dl2v1(1-p2)DlP∞(t)v1*p2DlP∞(t)*v1*v2Dl2(1-p1)p2(1-p1)v2DlP∞(t)v1*v2*Dl2[P∞(t)2]*v1*(1-p2)DlP∞(t)*(1-p1)v2*DlP∞(t)*(1-p1)(1-p2)],
where
(2.23)Dl=Dl(t)=e-ϰl2Γl(t),P∞(t)=e-iϰc2S(t)[1-2p]N1-2η.
The density matrix (2.22) is of the product form:
(2.24)[p1v1Dl(t)P∞(t)v1*Dl(t)P∞(t)*1-p1]⊗[p2v2Dl(t)P∞(t)v2*Dl(t)P∞(t)*1-p2].
This shows that the N→∞ asymptotic dynamics is factorizable and cannot create entanglement at any time.


## 2.1.1. Variation ofN Dependence

For homogeneous initial conditions,pj=p for j=3,…,N, we have
(2.15)PN(t)=[peiϰc2S(t)+(1-p)e-iϰc2S(t)]N-2.
Unless p=0,1, |PN(t)| oscillates in t between its minimum value |2p-1|N-2 (when cos(2ϰc2S(t))=-1) and its maximum value 1 (when cos(2ϰc2S(t))=1) (We have |peiϰc2S(t)+(1-p)e-iϰc2S(t)|2=p2+2p(1-p)cos(2ϰc2S(t))+(1-p)2.) The width of the oscillations becomes very narrow with increasing N. In the limit of large N, PN(t) is zero for all times, except for the discrete set of t∈ℝ satisfying cos(2ϰc2S(t))=1, in which case |PN(t)|=1.This implies that, for largeN, all off-diagonal density matrix elements of ρt vanish with the exception of [ρt]23=ei(ω1-ω2)te-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 PN(t)=1.)In order to try to have a nontrivial dynamics for largeN, one may scale the collective conserving coupling constant as
(2.16)ϰc⟶ϰcNη,someη>0.
Then (2.15) becomes
(2.17)PN(t)=[peiϰc2S(t)/N2η+(1-p)e-iϰc2S(t)/N2η]N-2.
An expansion in large N yields
(2.18)PN(t)=e-iϰc2S(t)[1-2p]N1-2ηe-2ϰc4S2(t)N1-4η[p(1-p)+O(N-2η)].
Thus as N→∞,
(2.19)PN(t)⟶{e-iϰc2S(t)[1-2p]N1/2e-2ϰc4S2(t)p(1-p)N1-4η,0<η≤14(p≠0,1),e-iϰc2S(t)[1-2p]N1-2η,14<η.Remark 2.3.
(1) By replacing, in these limits,ϰc2 by 2ϰc2, we obtain the corresponding limits for P~N(t).
(2) For0<η<1/4,  PN(t) vanishes as N→∞.
(3) The rapid oscillating phases disappear ifp=1/2 (any η>0) or η=1/2 (any p).


## 2.1.2. Asymptotic Concurrence (N→∞)

(
i
)
0
<
η
<
1
/
4
The reduced two-spin density matrix (in the interaction picture) at timet is
(2.20)[ρt′]=[p1p20000p1(1-p2)v1v2*e-2ϰl2Γl(t)00v1*v2e-2ϰl2Γl(t)(1-p1)p20000(1-p1)(1-p2)],
from which we obtain the concurrence
(2.21)C(ρt)=max⁡{0,-2[p1(1-p1)p2(1-p2)-|v1||v2|e-2ϰl2Γl(t)]}=0.
Remember that |vj|2≤pj(1-pj). This shows that the N→∞ asymptotic dynamics cannot create entanglement at any time.(
ii
)
η
>
1
/
4
Call the r.h.s. of (2.19) P∞(t) (a quantity still depending on N unless η=1/2). By replacing ϰc2 by 2ϰc2 in (2.19), 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→∞,
(2.22)[ρt′]=[p1p2p1v2DlP∞(t)v1p2DlP∞(t)v1v2Dl2P∞(t)2p1v2*DlP∞(t)*p1(1-p2)v1v2*Dl2v1(1-p2)DlP∞(t)v1*p2DlP∞(t)*v1*v2Dl2(1-p1)p2(1-p1)v2DlP∞(t)v1*v2*Dl2[P∞(t)2]*v1*(1-p2)DlP∞(t)*(1-p1)v2*DlP∞(t)*(1-p1)(1-p2)],
where
(2.23)Dl=Dl(t)=e-ϰl2Γl(t),P∞(t)=e-iϰc2S(t)[1-2p]N1-2η.
The density matrix (2.22) is of the product form:
(2.24)[p1v1Dl(t)P∞(t)v1*Dl(t)P∞(t)*1-p1]⊗[p2v2Dl(t)P∞(t)v2*Dl(t)P∞(t)*1-p2].
This shows that the N→∞ asymptotic dynamics is factorizable and cannot create entanglement at any time.


## 3. Numerical Results

(
i
) Let us first consider the case of two spins only, N=2. In (2.12) 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.8), we examine the maximal concurrence, as a function of S and Γ, for arbitrary fixed values of p and v.We find that, for fixedp, 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 Figure 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.(
ii
) 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.11), we choose the form factor, fc(k)=|k|χ|k|≤2πνc, with the cut-off frequency equal to the thermal frequency, νc=νT=kBT/h, at room temperature, T=300 K.In Figure2 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→ϰc2νct; see Figure 2(a), where νc is the cut-off frequency. The second effect is a reduction of the maximal concurrence in a smooth way; see Figure 2(b). 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 Figure 2(b).In Figure3 we show that the creation of concurrence decreases with the number of spins. In (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, Cmax, reported in (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 forN>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 Figure 3(c) and shows that the rescaled collapse time decays exponentially with the number of spins: τc=ϰc2νctc≃exp(-αN), where α=0.0838±0.0002. (See dashed line in Figure 3(c).)(
iii
) One can also vary the initial conditions for the spins by choosing independent p1,2 and v1,2, while all other spins have the same value pj=1/2, j=3,…,N (their off-diagonals vj do not influence the dynamics at all).In order to simplify the problem, we also setp1=v1 and p2=v2 and consider the maximal concurrence as a function of two independent parameters, p1 and p2, only.An example of the 3D plot obtained is reported in Figure4: the maximal concurrence is realized for p1=p2=v1=v2=1/2, independently of N (in the picture N=40, but similar plots are obtained for other values of N).(
iv
) 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 Figure5, 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 withN 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-off frequency in the range(0.5,4).


## 4. Conclusion

We have analyzed the two-qubit (two effective spin) entanglement in anN-qubit open system, where individual qubits interact 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 [16–18] 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 [19] 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 [20–23]. Here, the effective spins are exciton or electron sites (pigments), and the thermal environment represents the vibrational protein modes.



---

*Source: 375182-2012-07-10.xml*