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

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–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–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–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 section 2, we describe the model. In section 3, we discuss which observables we need to measure to obtain Bell inequality violations. In section 4, we present the results of our calculations in modeling the quantum dynamics of the system in time as quantum entanglement vanishes.

The target DQDs labeled A and B are described using a two-state basis composed of states $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{\alpha^A_0}=\Ket{0_A}$, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{\alpha^A_1}=\Ket{1_A}$, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{\alpha^B_0}=\Ket{0_B}$, and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{\alpha^B_1}=\Ket{1_B}$ in which the electron is fully localized on the bottom (0) 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 figure 1. Each of the environmental DQDs, indexed by $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} k=1, 2, \dots N$, is described by a similar two-state basis $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{\alpha^k_0}=\Ket{0_k}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} m\in \{0, 1\}$. Target DQD basis states are indexed with m_A and m_B, while the kth environmental DQD is indexed with m_k.

Zoom In Zoom Out Reset image size

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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \vec{m}\equiv [m_1, m_2, \dots, m_N]. \nonumber \end{align} \tag{ 1 }$$

There are $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} N_E=2^N$ such vectors, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} [\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 Ω, which includes the two target DQDs and the environment, consist of the direct products of the individual DQD states

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{\Phi_{m_A,m_B,\vec{m}_p}}= \quad \Ket{\alpha^A_{m_A} } \Ket{\alpha^B_{m_B} } \Ket{\alpha^1_{[\vec{m}_p]_1} } \Ket{\alpha^2_{[\vec{m}_p]_2} } \dots \Ket{\alpha^N_{[\vec{m}_p]_N} }. \label{eq:BasisPhis} \nonumber \end{align} \tag{ 2 }$$

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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{\psi(t)}= \sum_{{\begin{array}{@{}c@{}} m_A,m_B=0,1 \nonumber \\ 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 } }. \label{eq:TotalPsi} \nonumber \end{align} \tag{ 3 }$$

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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} U^{\,j, k}_{m_j, m_k}$ be the electrostatic potential energy between the jth qubit in state m_j (0 or 1) and the kth qubit in state m_k (0 or 1). This energy is given by

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle U^{\,j,k}_{m_j, m_k} = \frac{P(m_j) P(m_k) e^2}{16 \pi \epsilon_o} \left[ \frac{1}{r^{\,j,k}_{0,0}} - \frac{1}{r^{\,j,k}_{0,1}} - \frac{1}{r^{\,j,k}_{1,0}} + \frac{1}{r^{\,j,k}_{1,1}} \right]. \nonumber \end{align} \tag{ 4 }$$

where e is the fundamental charge, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \epsilon_o$ is the permittivity of free space, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} 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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} P(m)$ is the polarization of a DQD in state m. $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} P(1)=+1$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} 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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \vec{m}_p$ is calculated by simply summing over all interactions between all pairs of DQDs.

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle E_{m_A,m_B,\vec{m}_p} = \frac{1}{2}\sum_{j\ne k} U^{\,j,k}_{m_i,m_j} \; . \nonumber \end{align} \tag{ 5 }$$

Here the sums over indices i and j are over the DQDs $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} [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 equation (2) is:

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{H}=\sum_{m_A, m_B, p} \Ket{\Phi_{m_A,m_B,\vec{m}_p}} E_{m_A,m_B,\vec{m}_p} \Bra{\Phi_{m_A,m_B,\vec{m}_p}}. \label{eq:Htotal} \nonumber \end{align} \tag{ 6 }$$

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 0 to 1 with the environment in state $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \vec{m}_p$ as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle E^{{\rm {flip}}}_{A,\vec{m}_p}&\equiv E_{1,m_B,\vec{m}_p}-E_{0,m_B,\vec{m}_p},\nonumber \\ E^{{\rm {flip}}}_{B,\vec{m}_p}&\equiv E_{m_a,1,\vec{m}_p}-E_{m_A,0,\vec{m}_p}. \label{eq:FlipDef} \nonumber \end{align} \tag{ 7 }$$

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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} E^{{\rm {flip}}}_{{\rm {RMS}}}(A/B)$. A characteristic time can then be defined for each of the separated systems and the system as a whole as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \tau_{A,E} = h/E^{{\rm {flip}}}_{{\rm {RMS}}}(A), \quad \tau_{B,E} = h/E^{{\rm {flip}}}_{{\rm {RMS}}}(B), \quad \tau_E=\sqrt{\tau_{A,E} \tau_{A,E}}. \label{eq:TauDef} \nonumber \end{align} \tag{ 8 }$$

The characteristic time $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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 equation (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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{\rho}_\Omega(t)=\Ket{\psi(t)}\Bra{\psi(t)}\, = \sum_{{\begin{array}{@{}c@{}} m_A, m'_A \nonumber \\ m_B, m'_B \nonumber \\ p, p'\end{array}} } {c_{m'_A,m'_B,\vec{m}_p'} c^*_{m_A,m_B,\vec{m}_p } \Ket{\Phi_{m_A,m_B,\vec{m}_p }} \Bra{\Phi_{m'_A,m'_B,\vec{m}_{p'} }} }. \label{eq:RhoDef} \nonumber \end{align} \tag{ 9 }$$

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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \newcommand{\Argum}[1]{{\! \left( #1 \right)}} \displaystyle \hat{\rho}_\Omega \Argum{t} = {\rm e}^{-{\rm i} \frac{\hat{H}}{\hbar} t} \, \hat{\rho}_\Omega (0) \, {\rm e}^{+{\rm i} \frac{\hat{H}}{\hbar} t} . \label{eq:Liouville} \nonumber \end{align} \tag{ 10 }$$

This time evolution is exact within the model and the global system described by $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{\psi^{AB}(0)}=\left[ \Ket{0}_A\Ket{0}_B+\Ket{1}_A\Ket{1}_B \right]/\sqrt{2}. \label{eq:BellState} \nonumber \end{align} \tag{ 11 }$$

The initial state of the k_th environmental DQD is an unpolarized state given as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \newcommand{\MyPar}[1]{{\left( #1 \right)}} \newcommand{\Argum}[1]{{\! \left( #1 \right)}} \displaystyle \Ket{\psi_k \Argum{0}} = {\rm e}^{{\rm i}\theta_k} \MyPar{\Ket{0_k} + {\rm e}^{{\rm i} \phi_k} \Ket{1_k}}/\sqrt{2},\label{eq:initial_env} \nonumber \end{align} \tag{ 12 }$$

where the phases $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \theta_k$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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.

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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{0}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{1}$. In this basis we define the rotation operator $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \hat{R}$ as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{R}(\theta)=\cos(\theta)\left[\, \Ket{1}\Bra{1} + \Ket{0}\Bra{0}\,\right]+ \sin(\theta)\left[\, \Ket{0}\Bra{1} - \Ket{1}\Bra{0}\,\right]. \nonumber \end{align} \tag{ 13 }$$

We define two basis sets, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{a}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{a}'$, for measurements on subsystem A as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{a_l}&= \hat{R}(\theta_a) \Ket{l}_A, \nonumber \\ \Ket{a'_l}&= \hat{R}(\theta_{a'}) \Ket{l}_A \nonumber \end{align} \tag{ 14 }$$

for $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} l=0, 1$ for indicating the two basis states. We also define two basis sets, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{b}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{b}'$, for measurements on subsystem B

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{b_l}&= \hat{R}(\theta_b) \Ket{l}_B, \nonumber \\ \Ket{b'_l}&= \hat{R}(\theta_{b'}) \Ket{l}_B, \nonumber \end{align} \tag{ 15 }$$

for $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} l=0, 1.$ For the maximum Bell violation we choose $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} [\theta_a, \theta_{a'}, \theta_b, \theta_{b'}]=[0^{\circ}, 45^{\circ}, 22.5^{\circ}, 67.5^{\circ}]$. We define projection operators for measuring the four combinations of 0 and 1 on the two parties, for measurements using the $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{a}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \boldsymbol{b}$ bases as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{P}_{kl }(a,b)= \Ket{a_k}\Bra{a_k} \otimes \Ket{b_l}\Bra{b_l}, \label{eq:PplusminusDef} \nonumber \end{align} \tag{ 16 }$$

for $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} k, l=0, 1.$ The CHSH correlation function encodes the $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{0}$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \Ket{1}$ states with a  −1 and  +1 respectively. The expectation value for the product of the measurements (±1) on A and B using these bases is then given by

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \newcommand{\Tr}{{\rm Tr}} \displaystyle \hat{P}_{\times}(a,b)&\equiv P_{00}(a,b) - P_{01}(a,b)-P_{10}(a,b) + P_{11}(a,b), \nonumber \\ E(a,b)&=\left\langle P_{\times} \right\rangle =\Tr\left( \hat{\rho} \hat{P}_{\times}(a,b) \right). \label{eq:EchshDef} \nonumber \end{align} \tag{ 17 }$$

Expressions analogous to equations (16) and (17) define similar quantities $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} E(a, b')$, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} E(a', b)$, and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} E(a', b')$ using the other choices of basis states. The CHSH correlation function is then defined to be

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle S_{{\rm {CHSH}}}=\left| E(a,b) -E(a,b') + E(a',b) + E(a',b') \right|. \label{eq:CHSHcorrelation} \nonumber \end{align} \tag{ 18 }$$

The assumption of local realism yields the Bell inequality

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle S_{{\rm {CHSH}}} \leqslant 2. \label{eq:CHSHinequality} \nonumber \end{align} \tag{ 19 }$$

For local values of $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} a, a', b,$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} b'$ which are distributed randomly and uniformly $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {CHSH}}} = \sqrt{2}$ holds, obeying the inequality. For the fully entangled Bell state of equation (11), by contrast, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {CHSH}}} = 2\sqrt{2}$, in violation of (19).

3.2. Brukner–Paunković–Rudolph–Vedral correlation function

We will now consider the Bell inequality derived by Brukner et al in [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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \Ket{u_k}= \hat{R}(\theta_k) \Ket{0}, \quad \quad \Ket{v_k}= \hat{R}(\theta_k) \Ket{1}, \nonumber \end{align} \tag{ 20 }$$

for $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} k=1, 2, 3.$ For the correlation function we choose $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} [\theta_1, \theta_2, \theta_3]=[0^{\circ}, 120^{\circ}, 240^{\circ}]$. 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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{P}_k^{(0)} = \Ket{u_k}\Bra{u_k}, \quad \quad \hat{P}_k^{(1)} = \Ket{v_k}\Bra{v_k}, \nonumber \end{align} \tag{ 21 }$$

where we used the subscript $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} k=1, 2, 3$ for the three measurement settings. Each of these operators has eigenvalues 0 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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{P}_{kl}^{\rm 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), \label{eq:ProjectionDefs} \nonumber \end{align} \tag{ 22 }$$

and for obtaining the opposite outcomes as

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle \hat{P}_{kl}^{\rm 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), \label{eq:ProjectionDefs2} \nonumber \end{align} \tag{ 23 }$$

for $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} k\ne l.$ The operator $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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 0 or both 1).

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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle S_{{\rm {BPRV}}}= \sum^3_{k}\left \langle \hat{P}_{kk}^{\rm same}(A,B) \right\rangle + \sum^3_{k,l\ne k} \left \langle \hat{P}_{kl}^{\rm opp}(A,B) \right\rangle , \label{eq:MerminCorrelation} \nonumber \end{align} \tag{ 24 }$$

where the subscript refers to the initials of the authors of [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

$$\begin{align} \newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \displaystyle S_{{\rm {BPRV}}}\leqslant 7 \label{eq:MerminBellInequality} \nonumber \end{align} \tag{ 25 }$$

for any probability distribution of the measurement outcomes [37]. By contrast, the fully entangled Bell state of equation (11) yields $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {BPRV}}}=7.5$ in clear violation of equation (25). We evaluate directly the time-dependent value $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {BPRV}}}(t)$ from the density matrix evolving in time under equation (10).

Starting with the initial state given by equations (11) and (12), we solve for the unitary evolution of the global density matrix using equation (10), and calculate the correlation functions $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {CHSH}}}(t)$ and $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {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¹⁰  =  1024 environmental electronic configurations.

In figure 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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} 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.

Zoom In Zoom Out Reset image size

The CHSH correlation function $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} S_{{\rm {CHSH}}}(t)$ for each of these 72 configurations shown on the upper part of figure 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.

Figure 2(b) shows the same 72 cases of the geometrically random environment as figure 2(a), but plotted on a time axis scaled by the characteristic time $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \tau_E$ as calculated from equation (8). The value of $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \tau_{{\rm {opt}}}=1.34$ $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \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.

Figure 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.

Zoom In Zoom Out Reset image size

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 $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} t/\tau_E\approx 1.21$, whereas for the CHSH correlation function it occurs at $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} 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.

Figure 3(b) shows the dynamics of the von Neumann entropy of the reduced density operator for the two target qubits, $\newcommand{\Bra}[1]{\langle#1|}\newcommand{\Ket}[1]{|#1\rangle} \newcommand{\Tr}{{\rm Tr}}\hat\rho^{AB}(t)=\Tr_{E}(\hat\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.

Figure 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.

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).