Creation of Two-Particle Entanglement in Open Macroscopic Quantum Systems

Abstract

We consider an open quantum system of 𝑁 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 𝑁−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, 𝑁. Our results demonstrate that the two-qubit entanglement generally decreases as 𝑁 increases. We show analytically that, in the limit 𝑁→∞, 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.

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 of 𝑁 not directly interacting spins 1/2 (qubits) placed in a constant effective magnetic field (oriented in the 𝑧-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 𝑁−2 qubits. This gives us the time-dependent reduced density matrix for two arbitrary qubits. In the 𝑧-representation it is represented by a time-dependent 4×4 matrix. It is important to notice that the matrix elements, [𝜌_𝑡]_𝑛,𝑚, 𝑛,𝑚=1,…,4, depend not only on the parameters of the thermal environments but also on the total number of qubits, 𝑁. We study numerically the concurrence 𝐶(𝑡) of the reduced two qubit density matrix and its dependence on the parameters of the system and on 𝑁. 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, 𝐶_max(𝑡), generally decreases as 𝑁 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 the 𝑁 spins is in a state 𝜌_𝑗, 𝑗=1,…,𝑁; 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 Proposition 2.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: 𝜌_𝑡12=𝜌₀₁₂e^i𝜔2𝑡e^{i𝜘2𝑐𝑆(𝑡)}e^{−𝜘2ℓΓℓ(𝑡)−𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡). (1.1) The other matrix elements have similar behavior. Each factor on the r.h.s. has an interpretation.(i)[𝜌₀]₁₂ is the initial condition of the matrix element in question. None of the other initial matrix elements are involved (energy-conserving coupling).(ii)e^i𝜔2𝑡 is the uncoupled dynamics (no interaction with environments);(iii)e^{i𝜘2𝑐𝑆(𝑡)} is a dephasing factor with a time-dependent phase 𝑆(𝑡)≤0 becoming linear for large 𝑡 (for the considered infrared behavior |𝑘|^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 𝑁−2 traced-out spins (the term would be the same if only two spins were coupled to the reservoirs).(iv)e^{−𝜘2ℓΓℓ(𝑡)−𝜘2𝑐Γ𝑐(𝑡)} is a decaying factor with time-dependent decay rates, Γ(𝑡)≥0, becoming linear for large 𝑡 (see (2.11)). Both the local and collective reservoirs contribute; however, the term is independent of the 𝑁−2 traced-out spins (again, it would be the same if only two spins were coupled to the reservoirs);(v)𝑃_𝑁(𝑡) is a product of 𝑁−2 oscillating terms encoding the effect of all the traced-out spins (see (2.9)). It is important to notice that 𝑃_𝑁(𝑡) only depends on the diagonal density matrix elements of the initial states of the 𝑁−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 𝑁−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 𝑃_𝑁(𝑡) can be explained easily for the case in which all 𝑁−2 spins are initially in the high-temperature equilibrium state (1/√2)(|+⟩+|−⟩). Then 𝑃_𝑁(𝑡)=[cos(𝜘²_𝑐𝑆(𝑡))]^𝑁−2 and its magnitude oscillate between zero and one. Due to the large power (𝑁−2), the peaks of the function |𝑃_𝑁(𝑡)|, centered at the discrete times 𝑡_𝑝 satisfying 𝑆(𝑡_𝑝)∈(𝜋/𝜘²_𝑐)ℤ, are of very narrow width 𝑂(1/(𝜘²_𝑐√𝑁)) for large 𝑁. Consequently, in the limit 𝑁→∞, with 𝜘_𝑐 held fixed, |𝑃_𝑁(𝑡)| is zero for all 𝑡 except for 𝑡=𝑡_𝑝, where |𝑃_𝑁(𝑡_𝑝)|=1. But the density matrix becomes very simple if 𝑃_𝑁(𝑡)=0, because many entries vanish (c.f. Proposition 2.1) and the corresponding concurrence is zero. It follows that, in the large 𝑁 limit, concurrence is zero for all times (except possibly for some isolated instances, 𝑡_𝑝).

(ii) 𝑁-Dependent Scaling of the Interaction 𝜘_𝑐
The above analysis suggests that one cannot generate two-spin entanglement for 𝑁 large at fixed interaction strength 𝜘_𝑐. However, the width of the peaked function 𝑃_𝑁(𝑡) which is of order 𝑂(1/(𝜘²_𝑐√𝑁)) becomes appreciable if 𝜘_𝑐≳𝑁^1/4. Hence we consider an 𝑁-dependent scaling of the coupling, replacing 𝜘_𝑐 by 𝜘_𝑐/𝑁^𝜂, 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 𝜘_𝑐/𝑁^𝜂, we calculate the limit 𝑁→∞ of 𝜌_𝑡, for 𝑡∈ℝ 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 𝑡. While we are able to obtain explicit expressions for concurrence in the regime of 𝑁→∞, we are not so for 𝑁 finite. (The reduced density matrix is given explicitly for all 𝑁 and all 𝑡, but calculating from it explicitly the concurrence is more difficult.) However, since no entanglement is generated in the limiting case, 𝑁→∞, but we know entanglement is created for 𝑁=2 (see, e.g., [4, 11]), we expect that entanglement creation decays with increasing 𝑁. We study this decay numerically.

1.1.2. Numerical Results

We introduce 𝜈_𝑐, the highest frequency at which spin-reservoir interactions occur and call it the cut-off frequency. In the simulations, we take 𝜈_𝑐 of the order of the thermal frequency 𝜈_𝑇=𝑘_𝐵𝑇/ℎ. In the infra-red regime, our coupling is proportional to |𝑘| (see after (2.11)).

(i) For 𝑁=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 𝑁−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 𝑁.

Figure 1 

Maximal concurrence as a function of 𝑝 and 𝑣, for fixed Γ(𝑡)=0 ans 𝑆(𝑡)=𝜋/2. Here, 𝑁=2 spins is considered.

(ii) For general 𝑁, entanglement evolves according to a rescaled time 𝑡↦𝜘²_𝑐𝜈_𝑐𝑡, see Figure 2. This figure shows that a reduction of 𝜘_𝑐 diminishes the created concurrence in a moderate way. For instance, decreasing 𝜘_𝑐 by a factor 10 only decreases concurrence by less than 1/3.

(a)

(b)

Figure 2 

(a) Concurrence as a function of the rescaled time, for different 𝜘𝑐 and fixed 𝜖=𝜈𝑐/𝜈𝑇=1, being respectively the cut-off and thermal frequency (𝜈𝑇=𝑘𝐵𝑇/ℎ with 𝑇=300K) and the form factor 𝑓𝑐(𝑘)=|𝑘|. Other data are 𝑝=0.5,𝑣=0.48, 𝑁=4. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the coupling strength, 𝜘𝑐, for different 𝑁 values as indicated in the legend. Other values are the same as in (a).

(iii) In Figure 3 we show that the maximum of created concurrence decays with increasing 𝑁. For intermediate values of 𝑁 (with the current parameters N~10–150) the decrease is exponential; for smaller and larger values of 𝑁, it is superexponential.

(a)

(b)

(c)

Figure 3 

(a) Concurrence as a function of the rescaled time, for fixed 𝜖=𝜈𝑐/𝜈𝑇=1, being respectively the cut-off and thermal frequency (𝜈𝑇=𝑘𝐵𝑇/ℎ with 𝑇=300𝐾) and the form factor, 𝑓𝑐(𝑘)=|𝑘|. Other data are 𝑝=0.5, 𝑣=0.48, 𝜘𝑐=0.05. (b) Plot of the concurrence at the peak (obtained from (a)) as a function of the number of spins 𝑁. (c) Plot of the collapse time, 𝜏𝑐, as a function of the number of spins, 𝑁. The dashed line is the best fit exp(−𝛼𝑁) with 𝛼=0.0838±0.0002.

Figure 4 

Maximal concurrence as a function of the independent parameters, 𝑝1=𝑣1 and 𝑝2=𝑣2. As one can see the maximal concurrence is realized at the external corner, that is 𝑝1=𝑝2=1/2. Here is 𝜖=1, 𝑁=40 and 𝑝=1/2 for all other spins.

(iv) In the same Figure 3(c), we study the dependence of the maximal time, 𝜏_𝑐, (before recurrence) at which the concurrence is not zero. We have found that this time decays exponentially in the number of spins, 𝑁, for sufficiently large 𝑁.

(v) Results on the rescaled model 𝜘_𝑐↦𝜘_𝑐/𝑁^𝜂 are shown in Figure 5. We find a decrease of maximal concurrence with increasing 𝑁 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 𝑁, for intermediate values of 𝑁 (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 𝜘_𝑐↦𝜘_𝑐/𝑁^𝜂 can compensate the decay of created concurrence for large 𝑁.

Figure 5 

Maximal concurrence as a function of the number of spins, for different power law scaling, as indicated in the legend. Here is 𝑓𝑐(𝑘)=|𝑘|, 𝑝=0.5, 𝑣=0.48, 𝜖=𝜈𝑐/𝜈𝑇=1, 𝑇=300 K, 𝜘𝑐=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 the 𝑁 noninteracting spins 1/2 coupled by energy conserving interactions to local and collective bosonic heat reservoirs is given by 𝐻=−ħ^𝑁_𝑛=1𝜔_𝑛𝑆^𝑧_𝑛+^𝑁_𝑛=1𝐻_R𝑛+𝐻_R+^𝑁_𝑛=1𝜘_𝑛𝑆^𝑧_𝑛⊗𝜙_𝑐𝑓_𝑐+^𝑁_𝑛=1𝜈_𝑛𝑆^𝑧_𝑛⊗𝜙_𝑛𝑓_𝑛. (2.1)(2.2)

Below we use dimensionless variables and parameters. To do so, we introduce a characteristic frequency, 𝜔₀, typically of the order of spin transition frequency. The total Hamiltonian, energies of spin states, and temperature are measured in units ħ𝜔₀. The frequencies of spins, 𝜔_𝑛>0, bosonic excitations, 𝜔(𝑘)=𝑐|^→𝑘| (where 𝑐 is the speed of light), the wave vectors of bosonic excitations are normalized by 𝜔₀/𝑐, and all constants of interactions are measured in units of 𝜔₀. The dimensionless time is defined as 𝜔₀𝑡.

In (2.1) and (2.2), 𝜔_𝑛>0 is the frequency of spin 𝑛, 𝑆^𝑧=12100−1, (2.3) and 𝑆^𝑧_𝑛 denotes the 𝑆^𝑧 of spin 𝑛. 𝐻_R is the Hamiltonian of the bosonic collective reservoir, 𝐻_R=_ℝ3|𝑘|𝑎^∗(𝑘)𝑎(𝑘)d³𝑘, (2.4) and 𝐻_R𝑛 is that same Hamiltonian for the 𝑛th individual reservoir. For a square-integrable form factor ℎ(𝑘), 𝑘∈ℝ³, 𝜙(ℎ) is given by 𝜙(ℎ)=1√2^ℝ3ℎ(𝑘)𝑎^∗(𝑘)+ℎ(𝑘)^∗𝑎(𝑘)d³𝑘. (2.5) The real numbers, 𝜘_𝑛 and 𝜈_𝑛, 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 𝜘_𝑛=𝜘_𝑐(),𝜈∀𝑛collective_𝑛=𝜈_ℓ∀𝑛(local). (2.6) We also take, for simplicity, all local form factors equal (𝑓_ℓ) and all collective ones also (𝑓_𝑐).

Fix any pair of spins, and (re-)label their frequencies by 𝜔₁ and 𝜔₂; see (2.1). We write the reduced density matrix, 𝜌_𝑡, of the two fixed spins as a 4×4 matrix [𝜌_𝑡]_𝑖𝑗 in the ordered energy basis Φ₁=𝜑₁⊗𝜑₁,Φ₂=𝜑₁⊗𝜑₋₁,Φ₃=𝜑₋₁⊗𝜑₁,Φ₄=𝜑₋₁⊗𝜑₋₁, (2.7)(Another equivalent notation is: Φ₁=|++⟩,Φ₂=|+−⟩,Φ₃=|−+⟩,Φ₄=|−−⟩.) where 𝑆^𝑧𝜑_±1=±(1/2)𝜑_±1. For instance, [𝜌_𝑡]_2,4=⟨Φ₂,𝜌_𝑡Φ₄⟩.

The initial state of the spins is a product state of the form 𝜌_𝑆1,0⊗⋯⊗𝜌_𝑆𝑁,0, where 𝜌_𝑆𝑗,0=𝑝_𝑗𝑣_𝑗𝑣^∗_𝑗1−𝑝_𝑗, (2.8) with 0≤𝑝_𝑗≤1 and |𝑣_𝑗|²≤𝑝_𝑗(1−𝑝_𝑗). The upper bound on the off-diagonal guarantees that the eigenvalues of 𝜌_𝑆𝑗,0 are nonnegative.

We introduce the quantities: 𝑃_𝑁(𝑡)=^𝑁_𝑗=3𝑝_𝑗e^{i𝜘2𝑐𝑆(𝑡)}+1−𝑝_𝑗e^{−i𝜘2𝑐𝑆(𝑡)},𝑆(𝑡)1=−2_ℝ³𝑓_𝑐(𝑘)²|𝑘|𝑡−sin(|𝑘|𝑡)|𝑘|²d³Γ𝑘,_ℓ,𝑐(𝑡)=_ℝ3𝑓_ℓ,𝑐(𝑘)²𝛽|𝑘|coth2sin²(|𝑘|𝑡/2)|𝑘|²d³𝑘. (2.9)(2.10)(2.11)

The integrals in (2.10) and (2.11) are made to converge introducing a suitable cut-off wavenumber, |𝑘_𝑐|, or cut-off frequency, 𝜈_𝑐=|𝑘_𝑐|/2𝜋. (Here we use dimensionless units.) For instance, for numerical simulations, we choose as form factor the function 𝑓_𝑐(𝑘)=|𝑘|𝜒_{|𝑘|≤|𝑘𝑐|}, where 𝜒_{|𝑘|≤|𝑘𝑐|}=1 if |𝑘|≤|𝑘_𝑐| and 𝜒_{|𝑘|≤|𝑘𝑐|}=0 otherwise.

We also define 𝑃_𝑁(𝑡) to be the same as 𝑃_𝑁(𝑡), but with 𝜘²_𝑐 replaced by 2𝜘²_𝑐. 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 𝜌_𝑡12=𝜌₀₁₂e^i𝜔2𝑡e^{i𝜘2𝑐𝑆(𝑡)}e^{−𝜘2ℓΓℓ(𝑡)−𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡),𝜌_𝑡13=𝜌₀₁₃e^i𝜔1𝑡e^{i𝜘2𝑐𝑆(𝑡)}e^{−𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡),𝜌_𝑡14=𝜌₀₁₄e^{i(𝜔1+𝜔2)𝑡}e^{−2𝜘2ℓΓℓ(𝑡)−4𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡),𝜌_𝑡23=𝜌₀₂₃e^{i(𝜔1−𝜔2)𝑡}e^{−2𝜘2ℓΓℓ(𝑡)},𝜌_𝑡24=𝜌₀₂₄e^i𝜔1𝑡e^{−i𝜘2c𝑆(𝑡)}e^{−𝜘2ℓΓℓ(𝑡)−𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡),𝜌_𝑡34=𝜌₀₃₄e^i𝜔2𝑡e^{−i𝜘2𝑐𝑆(𝑡)}e^{−𝜘2ℓΓℓ(𝑡)−𝜘2𝑐Γ𝑐(𝑡)}𝑃_𝑁(𝑡), (2.12) and the populations are constant, [𝜌_𝑡]_𝑗𝑗=[𝜌₀]_𝑗𝑗, for 𝑗=1,…,4 and 𝑡∈ℝ.

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 spins 3,…,𝑁 is contained entirely in the factors 𝑃_𝑁(𝑡) and 𝑃_𝑁(𝑡). They only depend on the initial populations 𝑝_𝑗, 𝑗=3,…,𝑁 (see (2.9)), but not on the off-diagonals, 𝑣_𝑗. This is explained by the fact that when tracing over a single spin, 𝑗≥3, we perform the operation Tr_spin𝑗𝑈𝜌_𝑆𝑗,0𝑉, where 𝑈,𝑉 are operators commuting with 𝑆^𝑧₂ (energy-conserving interactions only!). Clearly the latter trace only involves the diagonal of 𝜌_𝑆𝑗,0.
(2) The oscillatory phases, e^i𝜔𝑡, in (2.12) represent the free, uncoupled dynamics of the spins. Consider the modified two-spin density matrix: 𝜌^′_𝑡=e^{i𝑡(−𝜔1𝑆𝑧1−𝜔2𝑆𝑧2)}𝜌_𝑡e^{−i𝑡(−𝜔1𝑆𝑧1−𝜔2𝑆𝑧2)}, (2.13) (“interaction picture” dynamics of 𝜌_𝑡). Because 𝜌^′_𝑡 and 𝜌_𝑡 are related by conjugation of a unitary operator of the product form, e^{i𝑡(−𝜔1𝑆𝑧1)}⊗e^{i𝑡(−𝜔2𝑆𝑧2)}, the concurrences of 𝜌_𝑡 and 𝜌^′_𝑡 are the same. In other words, when examining concurrence of 𝜌_𝑡, we may use formulas (2.12) with 𝜔₁=𝜔₂=0.

2.1. Concurrence

Recall that the concurrence of the reduced density matrix, 𝜌_𝑡, is unchanged when we pass to the interaction picture 𝜌^′_𝑡 (see the remark explaining (2.13)). In the basis (2.7), the evolution of 𝜌^′_𝑡 (2.13) is given by (2.12) with 𝜔₁=𝜔₂=0 and where the initial condition is (𝜌^′₀=𝜌₀): 𝜌₀=𝑝₁𝑝₂𝑝₁𝑣₂𝑣₁𝑝₂𝑣₁𝑣₂𝑝₁𝑣^∗₂𝑝₁1−𝑝₂𝑣₁𝑣^∗₂𝑣₁1−𝑝₂𝑣^∗₁𝑝₂𝑣^∗₁𝑣₂1−𝑝₁𝑝₂1−𝑝₁𝑣₂𝑣^∗₁𝑣^∗₂𝑣^∗₁1−𝑝₂1−𝑝₁𝑣^∗₂1−𝑝₁1−𝑝₂. (2.14)

2.1.1. Variation of 𝑁 Dependence

For homogeneous initial conditions, 𝑝_𝑗=𝑝 for 𝑗=3,…,𝑁, we have 𝑃_𝑁(𝑡)=𝑝e^{i𝜘2𝑐𝑆(𝑡)}+1−𝑝e^{−i𝜘2𝑐𝑆(𝑡)𝑁−2}. (2.15) Unless 𝑝=0,1, |𝑃_𝑁(𝑡)| oscillates in 𝑡 between its minimum value |2𝑝−1|^𝑁−2 (when cos(2𝜘²_𝑐𝑆(𝑡))=−1) and its maximum value 1 (when cos(2𝜘²_𝑐𝑆(𝑡))=1) (We have |𝑝e^{i𝜘2𝑐𝑆(𝑡)}+(1−𝑝)e^{−i𝜘2𝑐𝑆(𝑡)}|²=𝑝²+2𝑝(1−𝑝)cos(2𝜘²_𝑐𝑆(𝑡))+(1−𝑝)².) The width of the oscillations becomes very narrow with increasing 𝑁. In the limit of large 𝑁, 𝑃_𝑁(𝑡) is zero for all times, except for the discrete set of 𝑡∈ℝ satisfying cos(2𝜘²_𝑐𝑆(𝑡))=1, in which case |𝑃_𝑁(𝑡)|=1.

This implies that, for large 𝑁, all off-diagonal density matrix elements of 𝜌_𝑡 vanish with the exception of [𝜌_𝑡]₂₃=e^{i(𝜔1−𝜔2)𝑡}e^{−2𝜘2ℓΓℓ(𝑡)} (and [𝜌_𝑡]₃₂, of course) for almost all values of 𝑡. This suppression of off-diagonals comes from the large number of particles and is mediated through the collective energy-conserving interaction. (For 𝜘_𝑐=0 we have 𝑃_𝑁(𝑡)=1.)

In order to try to have a nontrivial dynamics for large 𝑁, one may scale the collective conserving coupling constant as 𝜘_𝑐⟶𝜘_𝑐𝑁^𝜂,some𝜂>0. (2.16) Then (2.15) becomes 𝑃_𝑁(𝑡)=𝑝e^{i𝜘2𝑐𝑆(𝑡)/𝑁2𝜂}+1−𝑝e^{−i𝜘2𝑐𝑆(𝑡)/𝑁2𝜂𝑁−2}. (2.17) An expansion in large 𝑁 yields 𝑃_𝑁(𝑡)=e^{−i𝜘2𝑐𝑆(𝑡)[1−2𝑝]𝑁1−2𝜂}e^{−2𝜘4𝑐𝑆2(𝑡)𝑁1−4𝜂[𝑝(1−𝑝)+𝑂(𝑁−2𝜂)]}. (2.18) Thus as 𝑁→∞, 𝑃_𝑁(𝑡)⟶e^{−i𝜘2𝑐𝑆(𝑡)[1−2𝑝]𝑁1/2}e^{−2𝜘4𝑐𝑆2(𝑡)𝑝(1−𝑝)𝑁1−4𝜂,}10<𝜂≤4,e𝑝≠0,1_{−i𝜘2𝑐𝑆(𝑡)[1−2𝑝]𝑁1−2𝜂},14<𝜂. (2.19)

Remark 2.3. (1) By replacing, in these limits, 𝜘²_𝑐 by 2𝜘²_𝑐, we obtain the corresponding limits for 𝑃_𝑁(𝑡).
(2) For 0<𝜂<1/4, 𝑃_𝑁(𝑡) vanishes as 𝑁→∞.
(3) The rapid oscillating phases disappear if 𝑝=1/2 (any 𝜂>0) or 𝜂=1/2 (any 𝑝).

2.1.2. Asymptotic Concurrence (𝑁→∞)

(i)0<𝜂<1/4
The reduced two-spin density matrix (in the interaction picture) at time 𝑡 is 𝜌^′_𝑡=𝑝₁𝑝₂0000𝑝₁1−𝑝₂𝑣₁𝑣^∗₂e^{−2𝜘2ℓΓℓ(𝑡)}00𝑣^∗₁𝑣₂e^{−2𝜘2ℓΓℓ(𝑡)}1−𝑝₁𝑝₂00001−𝑝₁1−𝑝₂, (2.20) from which we obtain the concurrence 𝐶𝜌_𝑡=max0,−2𝑝₁1−𝑝₁𝑝₂1−𝑝₂−𝑣₁𝑣₂e^{−2𝜘2ℓΓℓ(𝑡)}=0. (2.21) Remember that |𝑣_𝑗|²≤𝑝_𝑗(1−𝑝_𝑗). This shows that the 𝑁→∞ asymptotic dynamics cannot create entanglement at any time.

(ii)𝜂>1/4
Call the r.h.s. of (2.19) 𝑃_∞(𝑡) (a quantity still depending on 𝑁 unless 𝜂=1/2). By replacing 𝜘²_𝑐 by 2𝜘²_𝑐 in (2.19), we obtain the limit of 𝑃_𝑁(𝑡), which we call 𝑃_∞(𝑡). For 𝜂>1/4 we have the relation 𝑃_∞(𝑡)=[𝑃_∞(𝑡)]². The reduced two-spin density matrix (in the interaction picture) at time 𝑡 is, for 𝑁→∞, 𝜌^′_𝑡=𝑝₁𝑝₂𝑝₁𝑣₂𝐷_ℓ𝑃_∞(𝑡)𝑣₁𝑝₂𝐷_ℓ𝑃_∞(𝑡)𝑣₁𝑣₂𝐷²_ℓ𝑃_∞(𝑡)²𝑝₁𝑣^∗₂𝐷_ℓ𝑃_∞(𝑡)^∗𝑝₁1−𝑝₂𝑣₁𝑣^∗₂𝐷²_ℓ𝑣₁1−𝑝₂𝐷_ℓ𝑃_∞(𝑡)𝑣^∗₁𝑝₂𝐷_ℓ𝑃_∞(𝑡)^∗𝑣^∗₁𝑣₂𝐷²_ℓ1−𝑝₁𝑝₂1−𝑝₁𝑣₂𝐷_ℓ𝑃_∞(𝑡)𝑣^∗₁𝑣^∗₂𝐷²_ℓ𝑃_∞(𝑡)^2∗𝑣^∗₁1−𝑝₂𝐷_ℓ𝑃_∞(𝑡)^∗1−𝑝₁𝑣^∗₂𝐷_ℓ𝑃_∞(𝑡)^∗1−𝑝₁1−𝑝₂, (2.22) where 𝐷_ℓ=𝐷_ℓ(𝑡)=e^{−𝜘2ℓΓℓ(𝑡)},𝑃_∞(𝑡)=e^{−i𝜘2𝑐𝑆(𝑡)[1−2𝑝]𝑁1−2𝜂}. (2.23) The density matrix (2.22) is of the product form: 𝑝₁𝑣₁𝐷_ℓ(𝑡)𝑃_∞(𝑡)𝑣^∗₁𝐷_ℓ(𝑡)𝑃_∞(𝑡)^∗1−𝑝₁⊗𝑝₂𝑣₂𝐷_ℓ(𝑡)𝑃_∞(𝑡)𝑣^∗₂𝐷_ℓ(𝑡)𝑃_∞(𝑡)^∗1−𝑝₂. (2.24) This shows that the 𝑁→∞ 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, 𝑁=2. In (2.12) we put for simplicity Γ_ℓ=Γ_𝑐=Γ and regard Γ and 𝑆 as two independent parameters. Taking both spins initially in the same state given by 𝑝, 𝑣, see (2.8), we examine the maximal concurrence, as a function of 𝑆 and Γ, for arbitrary fixed values of 𝑝 and 𝑣.

We find that, for fixed 𝑝, 𝑣, the maximal concurrence is given at 𝑆=𝜋/2,Γ=0. Having such values fixed and plotting the concurrence as a function of 𝑝,𝑣, the maximal concurrence is realized when 𝑝=𝑣=1/2; see Figure 1, where a plot of the concurrence as a function of 𝑝,𝑣 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 𝑁 spins. For concreteness we choose 𝑝=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 𝑣=0.48. (Then 𝑝 is close to 𝑣 which favors larger entanglement creation.) Recall that the dynamics is independent of the off-diagonals of the 𝑁−2 traced-out spins (i.e., we do not have to specify the 𝑣 of the 𝑁−2 traced-out spins).

As mentioned after (2.11), we choose the form factor, 𝑓_𝑐(𝑘)=|𝑘|𝜒_{|𝑘|≤2𝜋𝜈𝑐}, with the cut-off frequency equal to the thermal frequency, 𝜈_𝑐=𝜈_𝑇=𝑘_𝐵𝑇/ℎ, at room temperature, 𝑇=300 K.

In Figure 2 we investigate the effects of an increase in the coupling parameter, 𝜘_𝑐. The first effect is a time shift for the concurrence evolution, described by a scaling, 𝑡→𝜘²_𝑐𝜈_𝑐𝑡; see Figure 2(a), where 𝜈_𝑐 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, 𝜘_𝑐, is relatively small. For instance, changing 𝜘_𝑐 for 𝑁=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 𝑁 values; see Figure 2(b).

In Figure 3 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 𝑁=2,…,32. As one can see, the same time rescaling is also valid for 𝑁>2. Moreover, the maximum concurrence created, 𝐶_max, reported in (b), decreases exponentially in 𝑁 in the range 10<𝑁<150 and faster than exponentially outside this range. For larger 𝑁, the concurrence decays superexponentially in 𝑁. For 𝑁 exceeding 200, the concurrence becomes too small to be significant (of the order 10⁻⁴).

It is also interesting to note that the graph of concurrence shows collapses and revivals and that the revival times for 𝑁>2 are always less than the revival time for 𝑁=2. It is also interesting to consider how the collapse time, 𝜏_𝑐, defined as the first time at which concurrence drops abruptly to zero, depends on the number of spins, 𝑁. This study has been reported in Figure 3(c) and shows that the rescaled collapse time decays exponentially with the number of spins: 𝜏_𝑐=𝜘²_𝑐𝜈_𝑐𝑡_𝑐≃exp(−𝛼𝑁), 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 𝑝_1,2 and 𝑣_1,2, while all other spins have the same value 𝑝_𝑗=1/2, 𝑗=3,…,𝑁 (their off-diagonals 𝑣_𝑗 do not influence the dynamics at all).

In order to simplify the problem, we also set 𝑝₁=𝑣₁ and 𝑝₂=𝑣₂ and consider the maximal concurrence as a function of two independent parameters, 𝑝₁ and 𝑝₂, only.

An example of the 3D plot obtained is reported in Figure 4: the maximal concurrence is realized for 𝑝₁=𝑝₂=𝑣₁=𝑣₂=1/2, independently of 𝑁 (in the picture 𝑁=40, but similar plots are obtained for other values of 𝑁).

(iv) The numerical analysis of the rescaled model with 𝜘_𝑐 replaced by 𝜘_𝑐/𝑁^𝜂 shows that the concurrence is always a decreasing function of 𝑁 and that the maximum of the created concurrence is a universal function of the number of spins 𝑁, independent of 𝜂 for 𝜂>1/4.

Results are shown in Figure 5, where the dashed line is the best exponential fit exp(−𝑎𝑁), with 𝑎=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 𝑁 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 an 𝑁-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 𝑁-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 𝑁. 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.

Acknowledgment

This work Supported by NSERC under Discovery Grant 205247 (M. Merkli). This work was carried out under the auspices of the NNSA of the U. S. DOE at LANL under Contract No. DEAC52-06NA25396. Supported by the LDRD program at LANL (G. P. Berman).