Quantum entanglement and transport in non-equilibrium interacting double-dot setup: The curious role of degeneracy

TL;DR

This paper investigates how quantum entanglement relates to electron transport in a non-equilibrium double-dot system, revealing the critical role of degeneracy and interactions in controlling entanglement levels.

Contribution

It demonstrates the influence of degeneracy and interactions on quantum entanglement in double-dot systems using NEGF, numerics, and asymptotic analysis.

Findings

Entanglement sharply decreases when bias exceeds a critical value in non-degenerate dots.

Degenerate dots maintain high entanglement, reaching 1/2 asymptotically.

Applying large chemical potentials can maximize entanglement in degenerate, non-interacting dots.

Abstract

We study quantum entanglement and its relation to transport in a non-equilibrium interacting double dot system connected to electronic baths. The dynamical properties in the non-interacting regime are studied using an exact numerical approach whereas the steady state properties are obtained following the well-known non-equilibrium Green's function (NEGF) approach. By means of mutual information and concurrence we explore the connection between the quantum correlations in the system and the current flowing through the dots. It is observed that entanglement between the dots is heavily influenced by the degeneracy or the lack thereof, of the dot levels. In the non-degenerate case, the concurrence falls sharply when the applied bias crosses a certain critical value. In contrast when the dot energy levels are degenerate, the concurrence reaches a very high asymptotic value of 1/2. When interactions are switched on, the degeneracy is lifted, and once again concurrence falls to zero beyond a critical value of the applied bias. Lastly it is observed that the concurrence can be made to reach almost the value of 1.0 if the chemical potential in both baths is made very large (while keeping the sign the same) provided the dot levels are kept degenerate within the non-interacting limit. A combination of NEGF method, brute-force numerics and asymptotics are employed to corroborate our findings.

Equations

$$H$$

$$H_D={ϵ}_1{n}_1+{ϵ}_2{n}_2+U{n}_1{n}_2.$$

$$H_{\alpha }=\sum _k{ϵ}_{k\alpha }{n}_{k\alpha },\alpha =L,R,$$

$$H_{\alpha D}=\sum _{i=1,2}{t}_{k\alpha }^i{c}_{k\alpha }^{†}{d}_i+h.c.,\alpha =L,R.$$

$${\Gamma }_{\alpha }^{ij}(ϵ)=2\pi \sum _k{t}_{k\alpha }^i(t_{k\alpha }^j{)}^{\ast }\delta (ϵ-{ϵ}_{k\alpha }).$$

$$G_{0,ij}(t,t^{\prime })=-i⟨T{d}_i(t)d_j^{†}(t^{\prime })⟩,i=1,2.$$

$$\Big{[}i\frac{\partial}{\partial t}-\epsilon_{i}\Big{]}G_{0,ij}(t,t^{\prime})=\delta_{ij}\delta(t-t^{\prime})+\sum_{k\alpha=L,R}(t^{i}_{k,\alpha})^{*}G^{0}_{kj}(t,t^{\prime}),$$

$$G_{kj}^0(t,t^{\prime })=-i⟨T{c}_{k\alpha }(t)d_j^{†}(t^{\prime })⟩$$

$$\Big{[}i\frac{\partial}{\partial t}-\epsilon_{k\alpha}\Big{]}G^{0}_{kj}(t,t^{\prime})=\sum_{l=1,2}t^{l}_{k\alpha}G^{0}_{lj}(t,t^{\prime}).$$

$$\displaystyle\Big{[}i\frac{\partial}{\partial t}-\epsilon_{i}\Big{]}g_{ij}(t,t^{\prime})$$

$$\displaystyle\Big{[}i\frac{\partial}{\partial t}-\epsilon_{k\alpha}\Big{]}g^{\alpha}_{k}(t,t^{\prime})$$

$${\Sigma }_{ij}^{\alpha }(t-t^{\prime })=\sum _k{t}_{k\alpha }^i{g}_{k\alpha }(t-t^{\prime })t_{k\alpha }^j,\alpha =L,R,$$

$${\mathbf{G}}_0(t,t^{\prime })=\mathbf{g}(t,t^{\prime })+{\int }^{t}d{t}_1{\int }^{t}d{t}_2\mathbf{g}(t,t_1)\mathbf{\Sigma }(t_1,t_2){\mathbf{G}}_0(t_2,t^{\prime }).$$

$${\mathbf{G}}_0(\tau ,{\tau }^{\prime })=\mathbf{g}(\tau ,{\tau }^{\prime })+\int d{\tau }_1\int d{\tau }_2\mathbf{g}(\tau ,{\tau }_1)\mathbf{\Sigma }({\tau }_1,{\tau }_2){\mathbf{G}}_0({\tau }_2,{\tau }^{\prime }).$$

$${\mathbf{G}}_0^{r/a}(ϵ)={\mathbf{g}}^{r/a}(ϵ)+{\mathbf{g}}^{r/a}(ϵ){\mathbf{\Sigma }}^{r/a}(ϵ){\mathbf{G}}_0^{r/a}(ϵ)$$

$${\bf G}_{0}^{r/a}(\epsilon)=\Big{[}\epsilon{\bf I}-{\bf H}_{D}-{\bf\Sigma}^{r/a}(\epsilon)\Big{]}^{-1}$$

$${\mathbf{G}}_0^{</>}(ϵ)={\mathbf{G}}_0^r(ϵ){\mathbf{\Sigma }}^{</>}(ϵ){\mathbf{G}}_0^a(ϵ),$$

$${\mathbf{\Sigma }}_{\alpha }^{r/a}(ϵ)=\mp i\frac{{\Gamma }_{\alpha }(ϵ)}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),$$

$${\mathbf{\Sigma }}_{\alpha }^{<}(ϵ)=i{f}_{\alpha }(ϵ){\Gamma }_{\alpha }(ϵ)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\text{and}$$

$${\mathbf{\Sigma }}_{\alpha }^{>}(ϵ)=-i(1-f_{\alpha }(ϵ)){\Gamma }_{\alpha }(ϵ)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).$$

$$t_{k\alpha }=\sqrt{\frac{{\Gamma }_{\alpha }({ϵ}_k)\Delta ϵ}{2\pi }}=\sqrt{\frac{{\Gamma }_{\alpha }\Delta ϵ}{2\pi }}.$$

$$H=\sum _{k=1}^N{\zeta }_k{c}_k^{†}{c}_k+\sum _{k=1}^{N-2}{t}_k(c_k^{†}{c}_{N-1}+c_{N-1}^{†}{c}_k)+t_k(c_k^{†}{c}_N+c_N^{†}{c}_k),$$

$$H=\left(\begin{array}{cccccccccc}{\zeta }_1& 0& \dots & 0& 0& 0& \dots & \dots & t_1& t_1\\ 0& {\zeta }_2& 0& \dots & \dots & \dots & \dots & 0& t_2& t_2\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ ...& ...& \dots & {\zeta }_{N_L}& 0& ...& ...& ...& t_{N_L}& t_{N_L}\\ 0& 0& \dots & 0& {\zeta }_1& 0& \dots & \dots & t_1& t_1\\ 0& 0& \dots & 0& 0& {\zeta }_2& \dots & \dots & t_2& t_2\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& \dots & 0& 0& 0& \dots & {\zeta }_{N_R}& t_{N_R}& t_{N_R}\\ t_1& t_2& \dots & t_{N_L}& t_1& t_2& \dots & t_{N_R}& {ϵ}_1& 0\\ t_1& t_2& \dots & t_{N_L}& t_1& t_2& \dots & t_{N_R}& 0& {ϵ}_2\end{array}\right).$$

$$\rho (0)={\rho }_L(0)\otimes {\rho }_R(0)\otimes {\rho }_D(0),$$

$${\rho }_L(0)=\frac{\exp (-\beta (H_L-{\mu }_L{\sum }_{k\in L}{c}_k^{†}{c}_k))}{Z_L},$$

$${\rho }_R(0)=\frac{\exp (-\beta (H_R-{\mu }_R{\sum }_{k\in R}{c}_k^{†}{c}_k))}{Z_R}.$$

$${\rho }_D(0)={\rho }_{D_1}\otimes {\rho }_{D_2}$$

$${\rho }_{D_1}(0)=n_0{d}_1^{†}{d}_1+(1-n_0)d_1{d}_1^{†}$$

$${\rho }_{D_2}(0)=n_0{d}_2^{†}{d}_2+(1-n_0)d_2{d}_2^{†},$$

$$\displaystyle I_{L/R}(t)=-\Big{\langle}\frac{dN_{L/R}}{dt}\Big{\rangle},$$

$$I_{L/R}=\int_{-\infty}^{\infty}\,\frac{d\epsilon}{2\pi}\,{\rm Tr}\Big{[}{\bf G}^{<}(\epsilon){\bf\Sigma}_{L/R}^{>}(\epsilon)-{\bf G}^{>}(\epsilon){\bf\Sigma}_{L/R}^{<}(\epsilon)\Big{]}.$$

$$I$$

$$I=\frac{1}{4}\int_{-\infty}^{\infty}\,\frac{d\epsilon}{2\pi}\,{\rm Tr}\Big{[}{\bf\Gamma}(\epsilon){\bf A}(\epsilon)\Big{]}\big{(}f_{L}(\epsilon)-f_{R}(\epsilon)\big{)}\Big{]},$$

Full text

I Introduction

Although an enormous body of knowledge has been created within the field of mesoscopic systems Beenakker (1991); Meir et al. (1991); Rokhinson et al. (1999); Matveev et al. (1996); Sztenkiel and Świrkowicz (2007); You and Zheng (1999); Zimbovskaya (2008), an understanding of the role of quantum correlations in these systems is rather primitive. Recent evidence Sharma and Rabani (2015); Sable et al. (2018) suggests that transport across quantum dots is intimately connected to the quantum correlations that develop between various subsystems. The question of how general these connections are, is still a matter of investigation. Non-equilibrium properties of interacting systems are notoriously hard to study, and there are few precious results on quantum correlations in such systems. Previous work Sharma and Rabani (2015); Sable et al. (2018) pertains only to the non-interacting limit, and moreover only mutual information (Wilms et al., 2011; Wolf et al., 2008; De Tomasi et al., 2017; Amico et al., 2008) has been studied, which contains not only quantum but also classical correlations. An understanding of the relationship between purely quantum correlations and transport, remains a pressing open question, even within the non-interacting limit. The other pressing question is the role of interactions. The present work is an attempt to make progress along both of these directions.

Here, we show how concurrence (Wootters, 1998; Cho and Kim, 2017; Deng et al., 2004; Zanardi and Wang, 2002; Nehra et al., 2018; Wu and Segal, 2011), which is an excellent measure of entanglement in fully mixed states may be studied in the presence of interactions. Furthermore our work pertains to the non-equilibrium regime, where we attempt a comparison of the transport properties and the quantum correlations that develop in the system. While the full time-evolution is accessible in the non-interacting limit, we concentrate on the steady state properties of the interacting model. Our system consists of a spinless quantum double dot connected to left and right baths with inter-dot interactions. In the current work we have theoretically studied the entanglement properties of a nonequilibrium interacting double dot system. A natural division of subsystems for the purpose of entanglement is the two dots, and we compute concurrence between the dots and investigate how it is connected to the current flowing through the dots.

Entanglement is very hard to measure experimentally. Quantum mechanics tells us that the expectation values of observables have the form $Tr(\rho O)$, where $O$ is the operator corresponding to the relevant observable. However, von Neumann entropy is given by $-Tr(\rho\log(\rho))$. $\log(\rho)$ is not an observable, and in fact, is a quantity that is dependent on the state itself. Therefore there are very subtle conceptual difficulties associated with whether it is even meaningful to measure von Neumann entropy. Hence, indirect access to entanglement is of great interest. In this context, there have been works that have tried to access entanglement via other quantities like current (Sharma and Rabani, 2015; Sable et al., 2018), conductance Yoo et al. (2018), and quantum noise Klich and Levitov (2009); Sharma and Rabani (2015). A lot of work has been done previously to understand the transport through the interacting quantum dot sytems and phenomena like the Coulomb blockade effect (Averin and Likharev, 1986; Waugh et al., 1995) and Kondo effect (Kondo, 1964) have been understood both theoretically and experimentally (Cronenwett et al., 1998; Ding et al., 2019; Goldhaber-Gordon et al., 1998; Brotons-Gisbert et al., 2019). Furthermore, within the Kondo regime, the entanglement between the magnetic impurity and the baths has been studied via different measures of entanglement (Yoo et al., 2018; Chung et al., 2018; Shim et al., 2018; Yang and Feiguin, 2017; Bonazzola et al., 2017). The indirect detection of entanglement via the conductance has also been proposed (Yoo et al., 2018). Clever strategies for a frontal attack on this problem continue to be of great current importance Pan et al. (2019).

The non-equilibrium Green’s Function (NEGF) method (Schwinger, 1961; Meir and Wingreen, 1992a; Haug and Jauho, 2008; Landauer, 1970; Wang et al., 2014; Jauho et al., 1994; Rammer and Smith, 1986) has been one of the most powerful analytical tools to tackle non-equilibrium physics. It has been used widely and successfully in the domain of quantum dots to study steady state transport and mesoscopic properties, and we too adopt this method. This technique is also useful in the non-interacting limit where we have access to the properties of the system in a deeper way. Here we take recourse to the equation of motion (EOM) technique (Levy and Rabani, 2013; Do, 2014; Li et al., 2012) to calculate the Green’s Function for both non-interacting and interacting models. For the non-interacting model we have followed both NEGF and exact numerical techniques in order to test for their mutual agreement, thus strengthening the quality of our results. For the interacting model we have used the Hartree-Fock approximation for decoupling the higher order correlators in such a way that the interaction effect on one site depends upon the average electron number of the other site. This is a simple mean field way of considering the interaction and yet has been argued to be legitimate for a large range of the Coulomb interaction strength (Meir et al., 1991; Sun and Guo, 2002; Levy and Rabani, 2013; Dhar and Sen, 2006; Agarwal and Sen, 2006; Lamba and Joshi, 2000; You and Zheng, 1999). There are also some advanced methods such as numerical renormalization group (NRG), real time renormalization group (RTRG), functional renormalization group (fRG), hierarchical quantum master equation (HQME) etc. (Lindner et al., 2019; Yoo et al., 2018; Erpenbeck et al., 2018; Härtle et al., 2013; Zheng et al., 2009; Jin et al., 2008; Okamoto et al., 2016; Härtle and Millis, 2014), but they are rather complicated and could be useful in extracting finer structure.

The study of quantum correlations within the non-interacting limit is facilitated by special tricks which exploit Wick’s Theorem, even for non-equilibrium systems Sharma and Rabani (2015). This allows for a reduction of a $2^{L}\times 2^{L}$ reduced density matrix problem to an $L\times L$ correlation matrix problem  Peschel (2003) with $L$ being the number of sites in the subsystem of interest, thus providing access to large system sizes. Exploiting this trick, it is possible to investigate mutual information between the quantum dots and the baths. However in the presence of interactions, the correlation matrix approach no longer holds, making the computation of quantities that involves reduced density matrix of the baths, a notoriously hard problem.

In the present work, we focus on the entanglement (mutual information, concurrence) that develops between two quantum dots when they are driven out-of-equilibrium, and investigate how it is correlated with the current flowing through the dots. The key findings of our paper are as follows: In the non-interacting limit, current and mutual information between the dots and the bath, have similar behaviour as in the single dot problem Sharma and Rabani (2015). The physics is strongly dependent on whether or not the two dot energy levels are degenerate. In the presence of degeneracy, both concurrence and mutual information between the dots, have finite steady state values, in contrast to their zero steady state values in the absence of degeneracy. The non-interacting degenerate system allows for the possibility of generating very high entanglement between the dots. By cranking up the left and right leads to a large and identical chemical potential, it becomes possible to make the two dots to reach unit concurrence value. Turning on interactions results in a lifting of the degeneracy leading to the concurrence again attaining a zero value in the steady state.

The plan of the paper is as follows. We begin with an introduction of the model Hamiltonian in Sec. II. In Sec. III we first discuss the non-interacting limit of the model following the NEGF approach and an exact numerical technique and present analytical and numerical results for charge transport, mutual information and concurrence. In Sec. IV we extend our study to the interacting Hamiltonian. Finally in Sec. V we summarize our findings.

II Model Hamiltonian

We consider an open quantum system with two parallel quantum dots connected to both left and right baths. The baths consist of non-interacting spinless fermions. In the non-interacting limit, the coherent tunnelling term between the dots is absent whereas in the interacting limit the dots exchange energy via Coulomb interaction. The total Hamiltonian can be separated into different components:

[TABLE]

where $H_{D}$ represents the dot Hamiltonian

[TABLE]

Here $n_{i}=d_{i}^{\dagger}d_{i},i=1,2$ represents the number operator with $d_{i}^{\dagger}$ $(d_{i})$ being the electronic creation (annihilation) operator for the $i$-th dot with energy $\epsilon_{i}$. $U$ denotes the Coulomb interaction strength.

[TABLE]

represents the bath Hamiltonian modelled as an infinite collection of non-interacting fermions with momentum index $k$ and corresponding energy $\epsilon_{k}$. $n_{\alpha k}=c^{\dagger}_{k\alpha}c_{k\alpha}$ is the number operator with $c_{k\alpha}^{\dagger}(c_{k\alpha})$ being the creation (annihilation) operators for $k^{\text{th}}$ mode. Finally, the system-bath coupling Hamiltonian $H_{\alpha D}$ is given as

[TABLE]

Here $t^{i}_{k\alpha}$ represents the tunneling amplitude between the $i^{\text{th}}$ dot and the $k^{\text{th}}$ mode of bath. This coupling information is encoded in the spectral density of the baths as follows:

[TABLE]

Throughout the paper we set $\hbar=e=k_{B}=1$. The dots are driven out-of-equilibrium by maintaining the electronic baths at different chemical potentials $\mu_{\alpha}$. In the present work, we keep the two bath temperatures identical $T_{L}=T_{R}$, although in principle this is another knob that can yield a non-equilibrium scenario. We then follow the dynamical as well as steady state properites of various thermodynamic and quantum informatic observables. In what follows, we will first consider the non-interacting limit ($U=0$) and present analytical and numerical results followed by results for the interacting case ($U\neq 0$).

III Non-interacting dots ($U=0$)

We begin with the NEGF formalism and discuss briefly the method to compute various two-point correlation functions (Green’s functions) for the subsystem as well as for the baths. In terms of these Green’s functions different thermodynamic observables can be computed.

III.1 NEGF Formalism

NEGF is a powerful tool to study transport properties of out-of-equilibrium many-body quantum systems. Within the NEGF method, there are several well established approaches to compute Green’s functions, like the Keldysh diagrammatic method Agarwalla et al. (2016); Benedict (2009), and the equation-of-motion (EOM) technique Niu et al. (1999); Levy and Rabani (2013); Galperin et al. (2007). Here we adopt the EOM approach to compute the Green’s functions. The first task is to define the time-ordered Green’s function for the dots:

[TABLE]

Here the time dependence in the operators represents the Heisenberg picture with operators evolving with the full Hamilonian $H$. $T$ represents the time-ordered operator and the average $\langle...\rangle$ is computed with respect to the initial condition. Following the Heisenberg’s EOM, the EOM for the Green’s functions can be obtained. A cascade of such equations of motion can then be written down, introducing new Green’s functions at every next step, until a closure among the Green’s functions is attained to terminate the procedure. For non-interacting systems, such a closure is possible to obtain which then yields exact analytical results. In contrast, for interacting systems, the closure is attained by employing different approximation schemes. The EOM for Eq. 6 is given as

[TABLE]

where

[TABLE]

is the time-ordered version of a mixed Green’s function involving the dot and the lead. In a similar manner, a differential equation for this mixed Green’s function can be obtained:

[TABLE]

Introducing the bare Green’s functions for the subsystem (dots) and the leads:

[TABLE]

and defining the self-energy for the electronic baths

[TABLE]

we can write down a formal solution for Eq. 7 as

[TABLE]

Here ${\bf\Sigma}={\bf\Sigma}_{L}+{\bf\Sigma}_{R}$ is the total self-energy, additive in both left and right leads, and is associated with the transfer of electrons between the leads and the dots. The bold notation refers to a matrix in the subsystem space. In a general non-equilibrium setup, a similar equation as in Eq. 12 can be obtained for the contour-ordered Green’s function where instead of real times $t$ and $t^{\prime}$ one introduces contour time variables $\tau$ and $\tau^{\prime}$ that runs on a complex time plane (Haug and Jauho, 2008). We therefore write, for non-equilibrium systems,

[TABLE]

From this contour-ordered Green’s function, applying Langreth theorem (Haug and Jauho, 2008), one then has access to all other real time Green’s functions namely, retarded $({\bf G}^{r}_{0}(t,t^{\prime}))$, advanced $({\bf G}^{a}_{0}(t,t^{\prime}))$, lesser $({\bf G}^{<}_{0}(t,t^{\prime}))$ and greater $({\bf G}^{>}_{0}(t,t^{\prime}))$ components. In the steady state limit, the expression for these Green’s functions can be simplified by taking advantage of the time-translational symmetry. In this limit, the Green’s functions depend only on the relative time difference $|t-t^{\prime}|$ and performing a Fourier transformation one obtains for retarded and advanced components

[TABLE]

which can be rewritten as

[TABLE]

with $H_{D}={\rm diag}(\epsilon_{1},\epsilon_{2})$. The lesser and greater components follow the Keldysh equation

[TABLE]

where ${\bf\Sigma}^{r,a,<,>}={\bf\Sigma}_{L}^{r,a,<,>}+{\bf\Sigma}_{R}^{r,a,<,>}$ are different components of the total self-energy given by:

[TABLE]

In writing the expressions for ${\bf\Sigma}_{\alpha}^{r/a}(\epsilon)$ we have ignored the real part which is responsible for the re-normalization of the bare dot energies. We further assume identical coupling between the bath and the dots $\Gamma^{ij}_{\alpha}=\Gamma_{\alpha}$. Note that the non-trivial temperature and bias information is contained in the lesser and greater components of the Green’s functions, given in Eq. 16.

III.2 Exact Numerical Approach and Dynamical Protocol

We next consider an exact numerical approach (Sharma and Rabani, 2015) to simulate the non-interacting double quantum dot setup. This method gives direct access to the density matrix of the entire system (dots, baths). In this approach, we discretize the bath spectrum and consider a large but equispaced finite number of levels. The coupling between the baths and the dots are then fixed by integrating the spectral function (Eq. 5) in a small energy window $\Delta\epsilon$ around $\epsilon_{k}$ and further by imposing the wide-band limit we obtain

[TABLE]

The final discretized form of the Hamiltonian $H$ can then be written as

[TABLE]

where $N$ is the total number of levels including the baths and the dots, $N=N_{L}+N_{R}+2$ with $N_{L}$ and $N_{R}$ being the energy-levels in the left and right baths, respectively. The indices $N-1$ and $N$ are used for the two dots with $c_{N-1}^{\dagger}(c_{N-1})=d_{1}^{\dagger}(d_{1}),c_{N}^{\dagger}(c_{N})=d_{2}^{\dagger}(d_{2})$, and the corresponding energies are, $\zeta_{N-1}=\epsilon_{1}$ and $\zeta_{N}=\epsilon_{2}$. The explicit form of the Hamiltonian in the matrix form can be expressed as

[TABLE]

Note that, one can write the above Hamiltonian $H$ in a diagonal form by introducing new set of fermionic operators $a_{\beta}$, defined as $a_{\beta}=\sum_{i=1}^{N}\psi_{\beta}(i)c_{i}$, where $\psi_{\beta}(i)$ are coefficients of the eigenvectors of the Hamiltonian coupling matrix.

Next we choose a decoupled initial condition for the global density matrix, written as a tensor product of the density matrix of each part as

[TABLE]

where the left and the right fermionic baths are distributed according to the grand canonical distribution

[TABLE]

Here both left and right baths are maintained at a fixed inverse temperature $\beta=1/T$ and at different chemical potentials $\mu_{L}$ and $\mu_{R}$ respectively. The density matrix for the dots is represented as :

[TABLE]

where $n_{0}$ is the initial population of the dots. Given this initial setup, the global density matrix at any arbitrary time is then simulated following the unitary time evolution $\rho(t)=e^{-iHt}\rho(0)e^{iHt}$.

For our numerical simulations, we have considered the baths with $N_{L/R}=128$ levels and finite cut-off for the energy spectrum $-\epsilon_{c}\leq\epsilon\leq\epsilon_{c}$ with $\epsilon_{c}=20$. We have fixed the step size as $\Delta\epsilon=2\epsilon_{c}/({N_{L}}-1)$. These parameters provide convergence upto the desired accuracy. Using this scheme we provide results for the current, mutual information and concurrence.

Note that, in the non-interacting case, the Gaussian nature of the initial state allows the use of Wick’s theorem to translate the problem of calculating the density matrix to a simpler problem of calculating the correlation matrix: $C_{ij}=\langle c_{i}^{\dagger}c_{j}\rangle$ (Sharma and Rabani, 2015; Dhar et al., 2012). Different physical observables can then be computed following this correlation matrix approach.

III.3 Results for charge current, mutual information and concurrence

We will begin by defining three relevant observables for this setup, namely the charge current, mutual information and concurrence.

Charge current: We first investigate the charge current flowing out of one of the leads. The current is defined as the rate of change of occupation number in a particular bath

[TABLE]

with $N_{L/R}=\sum_{k}c^{\dagger}_{kL/R}c_{kL/R}$. A formal expression for the steady state charge current for an arbitrarily interacting quantum junction can be written down in terms of the Green’s functions and is given by the celebrated Meir-Wingreen formula  (Haug and Jauho, 2008; Meir and Wingreen, 1992b)

[TABLE]

Here ${\bf G}^{</>}(\epsilon)$ refers to the Green’s function for an interacting subsystem and ${\bf\Sigma}_{L/R}^{</>}(\epsilon)$ is the usual self-energy term involving the leads. Upon symmetrization, $I=(I_{L}-I_{R})/2$, the expression for current simplifies to

[TABLE]

For a symmetric junction i.e., ${\bf\Gamma}_{L}(\epsilon)={\bf\Gamma}_{R}(\epsilon)={\bf\Gamma}(\epsilon)/2$ the first term of the above expression vanishes and one obtains a formula for steady state current for a general interacting system as

[TABLE]

where {\bf A}(\epsilon)=i\big{(}{\bf G}^{r}(\epsilon)-{\bf G}^{a}(\epsilon)\big{)} is defined as the spectral function matrix for the subsystem.

We now present the numerical results for the charge current. We display the steady state results obtained following the NEGF method whereas for the charge current dynamics we follow the exact numerical approach. Note that, for the numerical approach, we compute the following expression for current which can be obtained using Eq. 26 and the Heisenberg equation of motion,

[TABLE]

These correlators can be calculated for arbitrary time following the correlation matrix approach.

In Fig. (2a) and (3a) the dynamics of the charge current is shown both for degenerate and non-degenerate cases. After a short transient dynamics, in each case, the current saturates to a steady state value and this saturation value depends on the bias difference applied across the dots. However, interestingly, in the large bias limit, the asymptotic value is indifferent for both degenerate and non-degenerate dots and saturates to a value $1/2$.

We next compare the steady state current as a function of the bias voltage in Fig. (4a) and Fig. (5a). A nice agreement between the exact numerical approach and NEGF results is obtained for both degenerate and non-degenerate dots. The slight mismatch occurs due to the finite discretization of the bath spectrum. Note that, in the non-degenerate case, one observes an additional plateau in the current as compared to the degenerate case. Due to the symmetric choice of the chemical potential i.e., $\mu_{L}=-\mu_{R}$, the plateau starts to appear at $\mu_{L}-\mu_{R}\approx 2\epsilon_{1}$ providing a single channel for the electron to flow and disappears when $\mu_{L}-\mu_{R}\approx 2\epsilon_{2}$

Mutual information: To understand the corresponding effects in the context of correlations in the system, we study the mutual information (Wilms et al., 2011; Wolf et al., 2008; De Tomasi et al., 2017; Amico et al., 2008) between the dots and the baths as well as between the two dots. The mutual information for a bipartite system $A$ and $B$ is defined as

[TABLE]

where $S_{A}=-\text{Tr}_{A}(\rho_{A}\mbox{ln}\rho_{A})$, $S_{B}=-\text{Tr}_{B}(\rho_{B}\mbox{ln}\rho_{B})$ and $S_{AB}=-\text{Tr}_{A+B}(\rho_{AB}\mbox{ln}\rho_{AB})$ corresponds to the von Neumann entropies of $A,\ B$ and the composite system $A+B$, respectively. For example, to compute the mutual information between the two dots, $A$ and $B$ should correspond to individual dots.

Typically, for a general interacting system, the computation of the bath Von Neumann entropy is a challenging problem. However, in the non-interacting limit, following the correlation matrix approach (Peschel, 2003, 2012), the von Neumann entropy $S_{G}$ of a subspace $G$ can be computed as and is given as (Sharma and Rabani, 2015)

[TABLE]

where $C_{\sigma}$ are the eigenvalues of the correlation matrix defined within subspace $G$ and $N_{G}$ are the total number of sites in that subspace. For example, for subsystem consisting of two quantum dots $N_{G}=2$ and C_{ij}={\rm Tr}\Big{[}\rho(t)d_{i}^{\dagger}d_{j}\Big{]},i,j=1,2 following the notations used in Eq. II with $U=0$.

In Fig. (2b) and Fig. (3b) we display the dynamics of mutual information $s_{d}(t)$ between the two dots for both degenerate and non-degenerate cases. As can be seen, similar to the charge current, after a transient regime a steady state is achieved for the mutual information. It starts with a zero value due to the choice of our initial condition (see Eq. 25). For the same dot energies and in the large bias limit, the mutual information saturates to a finite value $\approx 0.43145$ whereas it vanishes for different dot energies. This dynamical behaviour is very similar also for the concurrence between the dots, as shown in Fig. (2c) and Fig. (3c) and a detailed discussion of this follows ahead. The steady state results for the mutual information are displayed in Fig. (4b) and Fig. (5b). The signature of the plateau in the current for the non-degenerate case (Fig. (5a)) gets reflected in the mutual information by a peak.

We also investigate the steady state mutual information between the dots and the bath as a function of the bias voltage in (Fig. 6). Interestingly, in this case, the qualitative behaviour is the same as that of the current flowing through the dots. In the non-degenerate case a plateau appears followed by a saturation value $4\,\text{ln}2$, whereas in the degenerate case the mutual information directly saturates to a value $2\,\text{ln}2$. These saturation values for various observables in the asymptotic limit can be obtained analytically and is discussed in the subsection D.

Concurrence: Although the exact numerical method allows to calculate the mutual information between the dot system and the baths in the non-interacting regime, a generalization to the interacting case is not trivial. However the double-dot system lends itself naturally to study two-site entanglement between the dots. A particularly useful quantity to measure this type of entanglement is called concurrence  (Wootters, 1998; Cho and Kim, 2017).

For a number conserving Hamiltonian the reduced density matrix for two sites can be written as (Nehra et al., 2018)

[TABLE]

where, $u=\langle(1-n_{1})(1-n_{2})\rangle,w_{1}=\langle(1-n_{1})n_{2}\rangle,w_{2}=\langle n_{1}(1-n_{2})\rangle,v=\langle n_{1}n_{2}\rangle\ \text{and}\ z=\langle d_{2}^{\dagger}d_{1}\rangle$. These various correlators can be calculated using Wick’s theorem.

The concurrence is then given by

[TABLE]

In Fig. (4c) and Fig. (5c) we display the steady state concurrence between the dots as a function of bias voltage. In the non-degenerate limit (Fig. (5c)), at low bias, concurrence starts with a very small value and suddenly falls off sharply to zero. This sharp fall occurs when $\mu_{L}-\mu_{R}\approx 2\epsilon_{1}$. This transition value of the bias voltage is plotted against the lower dot energy level $\epsilon_{1}$ (inset in Fig. (5c)), setting $\epsilon_{2}=10$. A linear dependence is obtained following both NEGF and exact numerics. For small bias, neither of the dot channels are activated and therefore the two dots can be present in a mixed entangled state. In this regime, an increasing bias voltage results in an increasing entanglement. Further increase in the bias voltage $\mu_{L}-\mu_{R}\geq 2\epsilon_{1}$ allows electron to tunnel through by populating the lowest energy level of the dot and keeping the higher energy dot empty. This leads to a separable state for the dots and concurrence drops to zero. This trend continues even when $\mu_{L}-\mu_{R}\geq 2\epsilon_{2}$ as the density matrix for the two dots always remain separable.

In contrast, for the degenerate limit, the concurrence value increases monotonically with applied bias Fig. (4c) before reaching an asymptotic steady state value of $1/2$ in the large bias limit. Because of degeneracy among the dot levels, an incoming electron from the bath can not differentiate between the two levels leading to a formation of a mixed entangled state. It further remains as an entangled state for higher bias voltage.

The other factors which can influence the concurrence are the temperature ($T$) and the coupling between the bath and the dots ($\Gamma$) as studied in Fig. 7. On varying the temperature, the concurrence also changes because of the thermal transport which initiates the transport of electrons even before the bias reaches the critical value. This also reflects in the concurrence as it can be seen that a finite correlation between the electrons of the two sites exists even for very low bias. The same argument holds on varying the coupling $\Gamma$. A greater value of the coupling means the transport is active even for small bias voltage and a large finite value of concurrence is obtained. It is also worth noting that the saturation value of concurrence is independent of both temperature and coupling.

III.4 Analytical results in the asymptotic limit

In this subsection we provide asymptotic large bias limit results for the charge current, mutual information, and concurrence.

We begin with the two point correlators $\langle d_{i}^{\dagger}d_{j}\rangle$ which are related to lesser or greater version of the corresponding Green’s function and can be computed by performing an integral in the energy domain

[TABLE]

In the asymptotic limit of large bias difference $(\mu_{L}\to\infty,\mu_{R}\to-\infty)$ the Fermi distribution reduces to $f_{L}(\epsilon)=1$ and $f_{R}(\epsilon)=0$. Further considering symmetric coupling and the wide-band limit i.e., ${\bf\Gamma}(\epsilon)={\bf\Gamma}$, the lesser component of the total self-energy (Eq. III.1) simplifies to,

[TABLE]

The Keldysh equation, given in Eq. 16, then reduces to

[TABLE]

One further receives, in this limit, the following relation using Eq. 15

[TABLE]

We can thus write

[TABLE]

The above integration over $A_{ii}$ yields $1$ and over $A_{ij}$ yields [math] for the non-degenerate case whereas for the degenerate case the integration over all the components gives the value $\frac{1}{2}$. Therefore, in the asymptotic limit, for the non-degenerate case

[TABLE]

and for the degenerate case,

[TABLE]

Charge current: The above analysis helps us to obtain the asymptotic expression for the steady state charge current for any interacting junction, given as

[TABLE]

both for degenerate and non-degenerate cases. This yields $I=1/2$ for $\Gamma=1$ which exactly matches with our numerical results as shown in Figs. ((2a) and (3a))

Mutual Information and Concurrence: We next analyze the asymptotic values for the mutual information and the concurrence. The reduced density matrix for the double-dot system is given as in Eq. 33. Using this, the reduced density matrix for the the subsystem, where dot $1$ is considered as one subsystem while the other dot as second subsystem, can be written as

[TABLE]

Calculating these traces we can write the reduced density matrices as

[TABLE]

The asymptotic values of different correlators in the non-degenerate case are (as calculated above)

[TABLE]

and for the degenerate case

[TABLE]

The asymptotic value of the mutual information between the two dots can then be calculated using Eq. 31, which vanishes for the non-degenerate case in this large bias limit, as also reflected in our numerics, Fig. 5b. In contrast, for the degenerate case the asymptotic value of mutual information is found to be $\approx 0.43145$, as also obtained numerically in Fig. 4b. In a similar manner, the asymptotic value of the mutual information between the bath and dot system can be computed. In the same large bias limit $(\mu_{L}\to\infty,\mu_{R}\to-\infty)$, and given that the initial dot occupancies are zero, i.e., $n_{1}=n_{2}=0$, (see Eq. 25), the overall system initially is in a pure state and therefore $S_{AB}=0$ and $S_{A}=S_{B}$, which imply that the mutual information between the bath and the system is simply related by $S=2S_{A}$. Now using Eq. 33, in the non-degenerate scenario one receives the asymptotic value as $4\,\text{ln}2$, whereas in the degenerate case it comes out to be $2\,\text{ln}2$ and matches exactly with Fig. 6.

The calculation of concurrence using Eq. 34 and Eqs. (45,46) is also straightforward. In the non-degenerate case since $|z|=0$, the concurrence vanishes whereas for the degenerate case it reaches an asymptotic value $1/2$ and further validates our numerical results as shown in Fig. 4c and Fig. 5c.

IV Interacting Dots ($U\neq 0$)

IV.1 Formalism

In this section we extend our study for the interacting case. The Hamiltonian for the spinless double quantum dot with inter-site Coulomb interaction is given in Eq.II. The NEGF calculation for the interacting model follows along the similar lines as the non-interacting case, however, in this case higher order Green’s functions needs to be computed and a suitable approximation scheme is therefore required to achieve closure in the EOM approach.

Let us first start writing the EOM for the interacting Green’s function for the dots,

[TABLE]

Where due to interaction higher order (two-particle) Green’s function appears, such as,

[TABLE]

and we need to compute the corresponding EOM as well,

[TABLE]

At this point, let us introduce the Green’s function associated only with the dots as,

[TABLE]

In the above EOM we use the steady state property, i.e., $\dot{n}_{m}(t)=0$ and $\langle n_{m}(t)\rangle$ is the average dot occupancy. This procedure yields a higher order correlator, mixing the dots and the leads

[TABLE]

So, a further EOM is necessary:

[TABLE]

In the above equation the Green’s functions (in the right hand side of the equation) for which $``l=m"$ vanish because multiple occupancy on a dot at the same time is forbidden due to Pauli exclusion principle. The mixed correlators containing operators from the lead and the dots such as $\langle c_{k\alpha}d_{i}^{{}^{\dagger}}\rangle$ are set to zero thus treating the coupling of the leads and system upto second order in $t_{k\alpha}$. Also the Green’s functions like Eq. 48 are decoupled using the mean field approach as $G_{{i}{m}{j}}(t,t^{\prime})=\langle n_{m}\rangle G_{ij}(t,t^{\prime})$. These two approximations have been shown to yield reliable results for a sufficiently high temperature (Meir et al., 1991; Sun and Guo, 2002; Levy and Rabani, 2013). This approximation leads to a closure for the cascade of Green functions in the EOM approach. As done for the non-interacting case, one can then write down an EOM for the contour-ordered version and obtain all the components in real time using Langreth’s thorem. The various correlators involving the system operators can then be computed self consistently using

[TABLE]

IV.2 Results

To understand the effect of Coulomb interaction on entanglement properties we first check the transport properties by computing the current and the corresponding differential conductance. We follow the NEGF method and use Eq. 29 to evaluate the current. In Fig. 8(a)) the standard Coulomb blockade effect is evident and the corresponding peaks are reflected in the differential conductance. In Particular, the blockade starts when the bias voltage corresponds to $2\epsilon_{1},\ 2\epsilon_{2},\ 2(\epsilon_{1}+U)$ and $2(\epsilon_{2}+U)$. Note that, once again the asymptotic value for the current is $1/2$.

Next, to compute the concurrence following Eq. 34 one needs to calculate different two-point and four-point correlators. Since for the interacting case correlators like $\langle n_{i}n_{j}\rangle$ can not be computed following the same trick as done in the non-interacting case, we employ here the Hartree-Fock approximation scheme  (Dhar and Sen, 2006; Agarwal and Sen, 2006) i.e.,

[TABLE]

Fig. (8b) and (8c) display the steady state concurrence for non-degenerate and degenerate case in presence of the interaction. For the non-degenerate case, a similar argument as the non-interacting model is admissible. For low bias, the concurrence is small and drops sharply to zero for $\mu_{L}-\mu_{R}>2\epsilon_{1}$. As evident from the numerics, the interaction has practically no effect on concurrence, in this case.

In contrast, in the degenerate scenario with finite interaction $U$ (Fig. 8c), the behaviour of concurrence is rich and can be analyzed in three distinct regimes for the bias voltage: (i) when the bias voltage is below the energy level of the dot i.e., $\mu_{L}-\mu_{R}\leq 2\epsilon$, the concurrence increases monotonically just like in the non-interacting case. (ii) For bias voltage in the range $2\epsilon\leq\mu_{L}-\mu_{R}\leq 2(\epsilon+U)$, the dot levels are non-degenerate due to the finite interaction $U$ and while one dot is occupied by electron the other electron needs to have an effective $2(\epsilon+U)$ energy to tunnel. Even if electron manages to tunnel through this level it leads to a separable state and entanglement between the two dots starts to decrease. (iii) further increase in the bias i.e., $\mu_{L}-\mu_{R}>2(\epsilon+U)$ the concurrence drops to zero as after filling one dot with energy $\epsilon$, the other electrons can tunnel through the shifted $2(\epsilon+U)$ level, leaving the two dots in a separable state. The inset of Fig. (8c) shows the plot for the critical bias, at which concurrence hits zero, with the energy level of the dot $\epsilon$ and a linear dependence is observed. In this figure the concurrence shows a step-like fall with the bias. This is due to the uneven change in the average electron number $\langle n_{i}\rangle$ with bias at higher values of bias, which in turn is a numerical artifact of the self consistently approximate computation of the average electron number.

In order to contrast transport and quantum correlations, we further investigate the situation when the two fermionic baths are maintained at the same chemical potential ($\mu_{L}=\mu_{R}=\mu$). Due to zero bias difference, net current into and out of the dots cancel out providing a dynamical zero, independent of the chemical potentials of the leads and the dot energy levels. However, under the same scenario, the entanglement between the dots depends on the value of the chemical potential, as well as the dot energies. In Fig. (9) we display the concurrence starting with the degenerate dot level case $(U=0)$. A high value for the concurrence is observed with increasing chemical potential in the non-interacting case and finally satures to the value 1.0 representing maximally entangled state for the dots. However, turning on the interaction lifts the degeneracy and for large chemical potential the two dots achieve a separable state leading to a zero concurrence. The inset in Fig. (9) shows the linear dependence of the chemical potential at which concurrence reaches zero with the dot energy level. Again the step-like fall in concurrence here is caused by the numerical artifacts in the self consistent calculation of average electron number.

V Summary

In summary, we have studied quantum transport and entanglement properties for out-of-equilibrium spinless parallel interacting quantum double dot setup. Employing NEGF and exact numerical approaches we have investigated quantum dynamics and the steady state properties for non-equilibrium charge current, mutual information and concurrence. It is found that both the transient and steady state behaviour of these observables is critically dependent on whether or not the dot energies are degenerate. In addition, strong correlations between these observables is found in both transient and steady state regimes. For example, In the non-degenerate case, for high bias, both the mutual information and concurrence approaches to zero in the steady state. Whereas for the degenerate case, a high value for both these observables is found to exist. In contrast, the asymptotic value for the current remain insensitive both in degenerate and non-degenerate limit. The appearance of plateau in the current for the non-degenerate case also reflected in mutual information between the dots as well as between the dots and the bath.

For the interacting case, we employ the Hartree-Fock approximation scheme to compute the steady state concurrence. The effects of interaction are once again tied up with degeneracy effects. The non-degenerate case is largely independent of interactions whereas in the degenerate case, the concurrence increases for small bias but drops to zero beyond a critical bias because of the lifting of the degeneracy. The characteristic value of this bias though, is dependent on the strength of interaction unlike in the non-degenerate case.

The future work will direct towards understanding the entanglement and transport properties for extended many-body quantum systems, in particular, scaling properties with the system size.

Acknowledgements

AD thanks Sreeraj Nair for helpful discussions. BKA thanks IISER-Pune for the start up funding. BKA and AS jointly acknowledge the International Centre for Theoretical Sciences for facilitating discussions during a visit for participating in the program - Indian Statistical Physics Community Meeting 2018 (Code: ICTS/Prog - ISPCM2018/02). AS is grateful to SERB which funds AD’s JRF via the startup grant (File Number: YSS/2015/001696). D.S.B acknowledges PhD fellowship support from UGC India.