Entanglement loss in molecular quantum-dot qubits due to interaction with the environment

TL;DR

This paper investigates how environmental interactions cause entanglement loss in molecular quantum-dot qubits, analyzing the dynamics of entanglement decay and its dependence on environmental geometry with implications for quantum computing.

Contribution

It models the entanglement dynamics of two quantum-dot qubits interacting with complex, randomly oriented environments, revealing universal scaling behaviors of entanglement loss.

Findings

Entanglement decays over time depending on environmental geometry.

Disentanglement times can be estimated with realistic molecular parameters.

Results suggest universal behavior in entanglement loss across different configurations.

Abstract

We study quantum entanglement loss due to environmental interaction in a condensed matter system with a complex geometry relevant to recent proposals for computing with single electrons at the nanoscale. We consider a system consisting of two qubits, each realized by an electron in a double quantum dot, which are initially in an entangled Bell state. The qubits are widely separated and each interacts with its own environment. The environment for each is modeled by surrounding double quantum dots placed at random positions with random orientations. We calculate the unitary evolution of the joint system and environment. The global state remains pure throughout. We examine the time dependence of the expectation value of the bipartite Clauser-Horne-Shimony-Holt (CHSH) and Brukner-Paunkovi\'c-Rudolph-Vedral (BPRV) Bell operators and explore the emergence of correlations consistent with local realism. Though the details of this transition depend on the specific environmental geometry, we show how the results can be mapped on to a universal behavior with appropriate scaling. We determine the relevant disentanglement times based on realistic physical parameters for molecular double-dots.

Figures

Equations

$$\vec{m}\equiv [m_1,m_2,\dots ,m_N].$$

$$|{\Phi }_{m_A,m_B,{\vec{m}}_p}⟩=|{\alpha }_{m_A}^A⟩|{\alpha }_{m_B}^B⟩|{\alpha }_{[{\vec{m}}_p{]}_1}^1⟩|{\alpha }_{[{\vec{m}}_p{]}_2}^2⟩\dots |{\alpha }_{[{\vec{m}}_p{]}_N}^N⟩.$$

$$\Ket{\psi(t)}=\sum_{\footnotesize{\begin{array}[]{c}m_{A},m_{B}=0,1\\ p=1,2,\dots N_{E}\end{array}}}c_{m_{A},m_{B},\vec{m}_{p}}(t)\Ket{\Phi_{m_{A},m_{B},\vec{m}_{p}}}.$$

$$U_{m_j,m_k}^{j,k}=\frac{P(m_j)P(m_k)e^2}{16\pi {ϵ}_o}\left[\frac{1}{r_{0,0}^{j,k}}-\frac{1}{r_{0,1}^{j,k}}-\frac{1}{r_{1,0}^{j,k}}+\frac{1}{r_{1,1}^{j,k}}\right].$$

$$E_{m_A,m_B,{\vec{m}}_p}=\frac{1}{2}\sum _{j\ne k}{U}_{m_i,m_j}^{j,k}.$$

$$\hat{H}=\sum _{m_A,m_B,p}|{\Phi }_{m_A,m_B,{\vec{m}}_p}⟩E_{m_A,m_B,{\vec{m}}_p}⟨{\Phi }_{m_A,m_B,{\vec{m}}_p}|.$$

$$E_{A,{\vec{m}}_p}^{\text{mbox}flip}$$

$$E_{B,{\vec{m}}_p}^{\text{mbox}flip}$$

$${\tau }_{A,E}=h/E_{\text{mbox}RMS}^{\text{mbox}flip}(A),{\tau }_{B,E}=h/E_{\text{mbox}RMS}^{\text{mbox}flip}(B),{\tau }_E=\sqrt{{\tau }_{A,E}{\tau }_{A,E}}.$$

$$\displaystyle\hat{\rho}_{\Omega}(t)=\Ket{\psi(t)}\Bra{\psi(t)}\,=\sum_{\footnotesize{\begin{array}[]{c}m_{A},m^{\prime}_{A}\\ m_{B},m^{\prime}_{B}\\ p,p^{\prime}\end{array}}}{c_{m^{\prime}_{A},m^{\prime}_{B},\vec{m}_{p}^{\prime}}c^{*}_{m_{A},m_{B},\vec{m}_{p}}\Ket{\Phi_{m_{A},m_{B},\vec{m}_{p}}}\Bra{\Phi_{m^{\prime}_{A},m^{\prime}_{B},\vec{m}_{p^{\prime}}}}}.$$

$${\hat{\rho }}_{\Omega }\left(t\right)=e^{-i\frac{\hat{H}}{\hslash }t}{\hat{\rho }}_{\Omega }(0)e^{+i\frac{\hat{H}}{\hslash }t}.$$

$$|{\psi }^{AB}(0)⟩=\left[{|0⟩}_A{|0⟩}_B+{|1⟩}_A{|1⟩}_B\right]/\sqrt{2}.$$

$$|{\psi }_k\left(0\right)⟩=e^{i{\theta }_k}\left(|0_k⟩+e^{i{\varphi }_k}|1_k⟩\right)/\sqrt{2},$$

$$\hat{R}(\theta )=\cos (\theta )\left[|1⟩⟨1|+|0⟩⟨0|\right]+\sin (\theta )\left[|0⟩⟨1|-|1⟩⟨0|\right].$$

$$|a_l⟩$$

$$|a_l^{\prime }⟩$$

$$|b_l⟩$$

$$|b_l^{\prime }⟩$$

$${\hat{P}}_{kl}(a,b)=|a_k⟩⟨a_k|\otimes |b_l⟩⟨b_l|,$$

$${\hat{P}}_{\times }(a,b)$$

$$E(a,b)$$

$$S_{\text{mbox}CHSH}=\left|E(a,b)-E(a,b^{\prime })+E(a^{\prime },b)+E(a^{\prime },b^{\prime })\right|.$$

$$S_{\text{mbox}CHSH}\le 2.$$

$$|u_k⟩=\hat{R}({\theta }_k)|0⟩,|v_k⟩=\hat{R}({\theta }_k)|1⟩,$$

$${\hat{P}}_k^{(0)}=|u_k⟩⟨u_k|,{\hat{P}}_k^{(1)}=|v_k⟩⟨v_k|,$$

$${\hat{P}}_{kl}^{\mathrm{same}}(A,B)={\hat{P}}_k^{(0)}(A)\otimes {\hat{P}}_l^{(0)}(B)+{\hat{P}}_k^{(1)}(A)\otimes {\hat{P}}_l^{(1)}(B),$$

$${\hat{P}}_{kl}^{\mathrm{opp}}(A,B)={\hat{P}}_k^{(0)}(A)\otimes {\hat{P}}_l^{(1)}(B)+{\hat{P}}_k^{(1)}(A)\otimes {\hat{P}}_l^{(0)}(B),$$

$$S_{\text{mbox}BPRV}=\sum _k^3⟨{\hat{P}}_{kk}^{\mathrm{same}}(A,B)⟩+\sum _{k,l\ne k}^3⟨{\hat{P}}_{kl}^{\mathrm{opp}}(A,B)⟩,$$

$$S_{\text{mbox}BPRV}\le 7$$

Full text

1 Introduction

Quantum entanglement plays a central role in quantum information science and quantum optics [1, 2]. Due to recent technological breakthroughs, it is now possible to create entangled states, for instance, with photons, cold trapped ions, ultracold atoms and solid state systems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Even if we successfully entangle two particles, they soon will become disentangled due to quantum correlations built up with the environment. Hence, environmentally induced loss of entanglement has received much recent attention.

Because the environment involves many degrees of freedom, most models adopt a phenomenological approach which includes, for example, coherence and energy decay times as model inputs. The dynamics of the entangled particles are then non-unitary [16, 17, 18, 19]. Few-particle systems can be more tractable yet still illuminate the behavior of much larger systems. For example, recent experiments have been successful at realizing a coherent system-environment dynamics with few particles, and saw thermal behavior emerge in a sub-system, even when the global system evolution is unitary with a fixed energy [20, 21].

Here, we consider a specific quantum system embedded in a few-particle environment for which we can calculate the global unitary system+environment dynamics exactly. The model environment is large enough for us to observe the disentangling dynamics of the embedded subsystem and see behavior similar to that expected from a large environment. In particular, we examine qubits realized by double quantum dots (DQD) with an extra electron, where the electron position encodes the quantum information. We study a system of two DQDs which are initially prepared in a maximally entangled state. These target DQDs are spatially separate and each interacts Coulombically with its own environment consisting of randomly placed similar DQDs with random orientations. Each environmental DQD interacts with the target DQD and also with all other DQDs in the same environment. Due to this interaction, the target DQDs become disentangled, which we follow by computing the Bell correlation function of Clauser-Horne-Shimony-Holt (CHSH) [22] and also that of Brukner-Paunković-Rudolph and Vedral (BPRV) [23], and look for the transition to classical (i.e., local realistic) behavior and the time scale on which this happens. Though the system is small enough to calculate the unitary dynamics exactly, it represents a large-dimensional Hilbert space, so aspects of the behavior of truly large systems emerge.

The intent of our model is not to propose a setup that can readily be realized in the laboratory. Rather we explore disentanglement in a concrete physically-motivated system that captures the key elements of the issue—separated but initially entangled pairs, and spatially distinct environments—in a tractable model amenable to direct numerical solution.

Recent experiments have entangled spatially separated double-dots using photons [24, 25], and phonons [26]. Remote electron spin systems have been entangled over a distance of more than a kilometer using microwave photons [27] and resulting Bell CHSH violations were measured. We do not concern ourselves here with the details of how the initial entanglement is established, but rather examine its decay due to entanglement with the environment.

The main characteristics of our model are the following.

(i) The joint state of system and environment remains pure during the dynamics. Thus, we model both the system and environment exactly. There are no stochastic or phenomenological terms added into the model.

(ii) The setup is physically motivated, rather than based on an abstract spin chain model. In fact, this system is a useful model for molecular mixed-valence double quantum dot systems and has the advantage of including a natural and physically realistic coupling mechanism [28, 29, 30, 31]. Such a double-dot system is a promising candidate for digital computing at the nanoscale [32].

(iii) The two target DQDs do not interact with each other and they have their own environments, which are separated and therefore not coupled to each other. This reflects the typical physical situation of spatially separated qubits and avoids artificial environment-mediated entanglement between the target systems.

(iv) There is Coloumbic coupling between the environment DQDs themselves, not only between the system DQDs and the environment DQDs. The coupling strength is computed from the distance and the orientation of the double dots.

(v) The model is the simplest possible to contain the necessary ingredients. The Hamiltonian consist of terms corresponding to the Coulomb energy only, while the electrons cannot tunnel between the dots. Hence the electrostatic energy remains constant during the dynamics. Moreover, in the computational basis, only the phases of the state vector components change, the amplitude remains constant.

(vi) Remarkably, even if a small number of environment DQDs are considered (we will show results for 10 environmental DQDs below), the decay of entanglement between the system DQDs results in a smooth decay of the Bell correlations. After appropriate normalization, all curves corresponding to various random arrangements of the environment double dots collapse to the same curve.

This model is an extension of a previous studies of the decoherence of a single double dot qubit state due to the environment [32, 33]. There, entanglement with the environmental drives the local system into Zurek “pointer states.” More complex internal dynamics have also been studied. Mixed valence molecules, which might realize $1$ nm size double-dot qubits, have additional nuclear motion to consider. Electron transfer from one dot to the other is coupled to vibrational modes of the nuclei, and ultimately the substrate [34]. In the present model we consider only rigid double-dots to focus on the issue of entanglement loss alone.

Our paper is organized as follows. In Sec. 2, we describe the model. In Sec. 3, we discuss which observables we need to measure to obtain Bell inequality violations. In Sec. 4, we present the results of our calculations in modeling the quantum dynamics of the system in time as quantum entanglement vanishes.

2 Model description

The target DQDs labeled $A$ and $B$ are described using a two-state basis composed of states $\Ket{\alpha^{A}_{0}}=\Ket{0_{A}}$, $\Ket{\alpha^{A}_{1}}=\Ket{1_{A}}$, $\Ket{\alpha^{B}_{0}}=\Ket{0_{B}}$, and $\Ket{\alpha^{B}_{1}}=\Ket{1_{B}}$ in which the the electron is fully localized on the bottom ([math]) or top ($1$) dot, respectively. A fixed charge of $+e/2$ resides at each dot, providing net charge neutrality for each DQD. The two initially entangled DQDs are far apart and each interact with a separate environment.

The environment is composed of several similar DQDs surrounding the target DQDs. The environmental DQDs are arranged in a sphere of radius $R$ with the positions on the sphere and the orientation of the double-dots chosen randomly as shown in Fig. 1. Each of the environmental DQDs, indexed by $k=1,2,\dots N$, is described by a similar two-state basis $\Ket{\alpha^{k}_{0}}=\Ket{0_{k}}$ and $\Ket{\alpha^{k}_{1}}=\Ket{1_{k}}$. There are $N/2$ environmental DQDs surrounding each target DQD ($N$ is always chosen to be even). We label individual DQD basis states with an integer $m\in\{0,1\}$. Target DQD basis states are indexed with $m_{A}$ and $m_{B}$, while the $k_{th}$ environmental DQD is indexed with $m_{k}$.

The electronic configuration associated with a specific environmental basis state can then be referred to using the vector

[TABLE]

There are $N_{E}=2^{N}$ such vectors, $[\vec{m}_{1},\vec{m}_{2},\vec{m}_{3},\dots\vec{m}_{N_{E}}]$, each representing a specific electronic configuration of the environment $E$.

The basis states for describing the global system $\Omega$, which includes the two target DQDs and the environment, consist of the direct products of the individual DQD states

[TABLE]

The global state can then be written as a linear combination of these basis states

[TABLE]

2.1 System Hamiltonian

The Hamiltonian for the global system, including $A$, $B$, and the environment $E$, is determined only by the electrostatic interaction between DQDs in the basis state electronic configurations. Let $U^{j,k}_{m_{j},m_{k}}$ be the electrostatic potential energy between the $j^{\mbox{\tiny{th}}}$ qubit in state $m_{j}$ (0 or 1) and the $k^{\mbox{\tiny{th}}}$ qubit in state $m_{k}$ (0 or 1). This energy is given by

[TABLE]

where $e$ is the fundamental charge, $\epsilon_{o}$ is the permittivity of free space, $r^{j,k}_{m_{j},m_{k}}$ is the distance between dot $m_{j}$ in DQD $j$ and dot $m_{k}$ in DQD $k$, and $P(m)$ is the polarization of a DQD in state $m$. $P(1)=+1$ and $P(0)=-1$.

The total electrostatic potential energy of a configuration of target DQDs in states $m_{A}$ and $m_{B}$, and the environment in the state defined by $\vec{m}_{p}$ is calculated by simply summing over all interactions between all pairs of DQDs.

[TABLE]

Here the sums over indices $i$ and $j$ are over the DQDs $[A,B,1,2,\dots,N]$, that is, including both target and environmental DQDs.

In this model there is no tunneling between dots within the DQD; we are interested in the entanglement of the phase degrees of freedom rather than electron transfer effects which have been studied elsewhere [34]. The Hamiltonian for the global system is then diagonal in the basis states defined by Eq. (2) is:

[TABLE]

We can characterize the strength of the interaction between each target DQD and its local environment by the electrostatic energy needed to flip DQD $A$ or $B$ from [math] to $1$ with the environment in state $\vec{m}_{p}$ as

[TABLE]

These energies depend on the electrostatic configurations of the environment DQDs. We now define a quantity independent of the quantum state of the environment. Let the root-mean-square of the flip energies over all the electronic configurations of the local environmental basis states be $E^{\mbox{\tiny{flip}}}_{\mbox{\tiny{RMS}}}(A/B)$. A characteristic time can then be defined for each of the separated systems and the system as a whole as

[TABLE]

The characteristic time $\tau_{E}$ depends on the details of the geometrically random orientation and positions of the environmental DQDs. As we will see later, while the system is described by a complicated interaction of randomly placed double dots, the time constants given in Eq. (8) characterize the main aspects of the dynamics [33].

2.2 Density operator and equation of motion

The time evolution of the system is calculated using the equation of motion of the density operator. The density operator for the global system is defined from (3) by

[TABLE]

The dynamics of the global system density matrix is obtained by solving the von Neumann equation as

[TABLE]

This time evolution is exact within the model and the global system described by $\hat{\rho}_{\Omega}$ is always in a pure state.

We now define the initial state of the joint system. The target DQDs $A$ and $B$ are initially in the symmetric entangled Bell state

[TABLE]

The initial state of the $k_{th}$ environmental DQD is an unpolarized state given as

[TABLE]

where the phases $\theta_{k}$ and $\phi_{k}$ are chosen randomly, with a distribution that results in the corresponding Bloch vectors being uniformly distributed over the unit sphere. We take the initial state of the density operator to be a tensor product state of the entangled system $AB$ and the complete environment.

3 Tracking entanglement with Bell operators

We will primarily observe the disentanglement of the target DQDs by computing the dynamics of the expectation values of Bell operators. These are are relevant experimentally, since they can obtained by projective measurements on the subsystems. In the next section we will also calculate entanglement measures for the evolving system.

3.1 CHSH correlation function

We now calculate the Clauser-Horne-Shimony-Holt (CHSH) correlation function and the corresponding Bell inequality [22]. This function is often measured experimentally and it has been shown that states violating the CHSH inequality can be used in the Ekert protocol for entanglement assisted quantum communication [35].

For each subsystem $A$ and $B$ we define operators in the space spanned by the local basis vectors $\Ket{0}$ and $\Ket{1}$. In this basis we define the rotation operator $\hat{R}$ as

[TABLE]

We define two basis sets, $\bm{a}$ and $\bm{a}^{\prime}$, for measurements on subsystem $A$ as

[TABLE]

for $l=0,1$ for indicating the two basis states. We also define two basis sets, $\bm{b}$ and $\bm{b}^{\prime}$, for measurements on subsystem $B$

[TABLE]

for $l=0,1.$ For the maximum Bell violation we choose $[\theta_{a},\theta_{a^{\prime}},\theta_{b},\theta_{b^{\prime}}]=[0$\mathrm{\SIUnitSymbolDegree}$,45$\mathrm{\SIUnitSymbolDegree}$,22.5$\mathrm{\SIUnitSymbolDegree}$,67.5$\mathrm{\SIUnitSymbolDegree}$]$. We define projection operators for measuring the four combinations of [math] and $1$ on the two parties, for measurements using the $\bm{a}$ and $\bm{b}$ bases as

[TABLE]

for $k,l=0,1.$ The CHSH correlation function encodes the $\Ket{0}$ and $\Ket{1}$ states with a $-1$ and $+1$ respectively. The expectation value for the product of the measurements ($\pm 1$) on $A$ and $B$ using these bases is then given by

[TABLE]

Expressions analogous to Eq.(19) and (20) define similar quantities $E(a,b^{\prime})$, $E(a^{\prime},b)$, and $E(a^{\prime},b^{\prime})$ using the other choices of basis states. The CHSH correlation function is then defined to be

[TABLE]

The assumption of local realism yields the Bell inequality

[TABLE]

For local values of $a,a^{\prime},b,$ and $b^{\prime}$ which are distributed randomly and uniformly $S_{\mbox{\tiny{CHSH}}}=\sqrt{2}$ holds, obeying the inequality. For the fully entangled Bell state of Eq. (14), by contrast, $S_{\mbox{\tiny{CHSH}}}=2\sqrt{2}$, in violation of (22).

3.2 Brukner-Paunković-Rudolph-Vedral correlation function

We will now consider the Bell inequality derived by Č. Brukner et al., in Ref. [23], which is a generalization of Mermin’s Bell inequality [36, 37] for bipartite states that are not necessarily symmetric.

We define three sets of rotated basis function as

[TABLE]

for $k=1,2,3.$ For the correlation function we choose $[\theta_{1},\theta_{2},\theta_{3}]=[0$\mathrm{\SIUnitSymbolDegree}$,120$\mathrm{\SIUnitSymbolDegree}$,240$\mathrm{\SIUnitSymbolDegree}$]$. For each subsystem $A$ and $B$ a particular basis set corresponding to one of these angles is randomly chosen and a projective measurement is carried out. The projection operators corresponding to measurements are

[TABLE]

where we used the subscript $k=1,2,3$ for the three measurement settings. Each of these operators has eigenvalues [math] and $1$, hence all measurements on either subsystem $A$ or $B$ have these outcomes. We define the correlation operators for obtaining the same outcome in the two qubits as

[TABLE]

and for obtaining the opposite outcomes as

[TABLE]

for $k\neq l.$ The operator $\hat{P}_{12}^{\rm same}(A,B)$, for example, corresponds to a measurement of setting $1$ on one subsystem and setting $2$ on the other, both yielding the same result (both [math] or both $1$).

The correlation function needed to evaluate the Bell inequality is then defined as

[TABLE]

where the subscript refers to the initials of the authors of Ref. [23]. Local realism requires that each subsystem $A$ and $B$ have values that determine the results of measurements of settings $1$, $2$, and $3$ before the measurement is made. That assumption yields the inequality

[TABLE]

for any probability distribution of the measurement outcomes [37]. By contrast, the fully entangled Bell state of Eq. (14) yields $S_{\mbox{\tiny{BPRV}}}=7.5$ in clear violation of Eq. (28). We evaluate directly the time-dependent value $S_{\mbox{\tiny{BPRV}}}(t)$ from the density matrix evolving in time under Eq. (13).

4 Results

Starting with the initial state given by Eqs. (14) and (15), we solve for the unitary evolution of the global density matrix using Eq. (13), and calculate the correlation functions $S_{\mbox{\tiny{CHSH}}}(t)$ and $S_{\mbox{\tiny{BPRV}}}(t)$ directly from the global density matrix. The number of environmental double-dots is $N_{E}=10$ (five around each target double-dot) yielding $2^{10}=1024$ environmental electronic configurations.

In Fig. 2(a), the results are shown for $a=1$ nm, a typical scale for molecular double-dots, which sets the time scale at picoseconds. Calculated Bell correlation functions are shown for 6 different values of $R/a\in\{2.5,3,3.5,4,5,7\}$, corresponding to different average strengths of coupling to the environment. For each value of $R/a$, 12 different random geometric arrangements of the environments are shown.

The CHSH correlation function $S_{\mbox{\tiny{CHSH}}}(t)$ for each of these 72 configurations shown on the upper part of Fig. 2, begins at the value corresponding to maximally violating the Bell inequality. Then, it starts dropping out of the Bell violation regime and decaying to the classical limit. As expected, the stronger the coupling to the environment, the faster the quantum entanglement disappears.

Fig. 2(b) shows the same 72 cases of the geometrically random environment as Fig. 2(a), but plotted on a time axis scaled by the characteristic time $\tau_{E}$ as calculated from Eq. 8. The value of $\tau_{E}$ is distinct for each of the random geometries of the environment. The time-scaled result is independent of the values of $a$ or $R/a$.

The squares show the value of a single Gaussian fit to all 72 curves for the transition from the initial Bell-state value to the classical limit. The fit yields a Gaussian width of $\tau_{\mbox{\tiny{opt}}}=1.34$ $\tau_{E}$ and matches the calculated bundle of trajectories well. The CHSH correlation function evolves from the fully entangled value to the value corresponding to local realism over a time on the scale of $\tau_{E}$ and the transition is very close to Gaussian, rather than the often-assumed exponential associated with semigroup behavior. Note that the slope at small times is here zero, in contrast to the finite slope of an exponential.

Fig. 3(a) shows the time-scaled BPRV correlation, which behaves similarly, crossing out of the Bell violation regime and into the classical (local realism) limit with a Gaussian form. The width of the Gaussian is identical to that for the CHSH correlation function; the squares on the plot show the Gaussian fit. It may be that as the number of environmental DQDs increases, the phase interference that results in the slight residual oscillations averages out to a precisely Gaussian shape.

Though both correlation functions have the same Gaussian shape and the same width, the transition out of the region forbidden by local realism occurs at different times for the two different Bell correlation functions. For the BPRV correlation function, the transition occurs at $t/\tau_{E}\approx 1.21$, whereas for the CHSH correlation function it occurs at $t/\tau_{E}\approx 1.78$. Of course entanglement can persist even after the system no longer violates a particular Bell inequality.

The Gaussian shape of the transition into local realism is notable. A similar Gaussian characteristic has been observed by Cucchietti et al. [38] in the context of the decoherence of a model spin system. This result seems to hold across many distributions of coupling to the environment and is rooted in the approximately Gaussian distribution of the eigenvalues of the Hamiltonian for a random environment.

Fig. 3(b) shows the dynamics of the von Neumann entropy of the reduced density operator for the two target qubits, $\rho^{AB}(t)=\Tr_{E}(\rho_{\Omega}(t))$. The increase in entropy by one bit corresponds to the loss of local information about the state of the pair. The global von Neumann entropy remains zero throughout because the global state is always pure.

Fig. 3 also shows the time dependence of the entanglement of formation for the two target qubits [39]. This quantity decreases as the entanglement between the two DQDs smoothly vanishes.

5 Conclusion

We have examined the dynamics of qubits realized with a double-dot with an extra electron. Initially, they were placed in a maximally entangled symmetric Bell state. Then, each interacted Coulombically with its own environment consisting of similar double-dots with a random position and orientation. While the dynamics was unitary and the entire system remains pure, the two system double-dots become gradually disentangled from each other. We followed this by calculating the dynamics of the expectation values of Bell operators for relevant Bell inequalities. We used the physical parameters of molecular mixed-valence double quantum dot systems for modeling, and calculated the relevant times scales. While the dynamics were different for different random geometries, with appropriate normalization all curves collapsed to a single curve. The time scale of disentanglement could be calculated precisely from the mean interaction strength between the system and the environmental degrees of freedom.

Our results can be understood by noting that quantum entanglement is best characterized not as fragile, but rather as promiscuous. Dynamics entangles each system with all the other systems with which it interacts. This promiscuity is constrained by the principle of the quantum monogamy of entanglement [40, 41], which bounds the strength of entanglement between any two pairs when a system entangles with many other systems. The decay of entanglement we see between the two target DQD systems occurs precisely because they each entangle with multiple systems in their respective environments, while maintaining global coherence and purity completely.

We acknowledge support from the National Science Foundation (Grant No. DGE-1313583). We also acknowledge the support of the EU (ERC Starting Grant 258647/GEDENTQOPT, COST Action CA15220, QuantERA CEBBEC), the Spanish Ministry of Economy, Industry and Competitiveness and the European Regional Development Fund FEDER through Grant No. FIS2015-67161-P (MINECO/FEDER, EU), the Basque Government (Grant No. IT986-16), and the National Research, Development and Innovation Office NKFIH (Grant No. K124351).