Open system approach to non-equilibrium dynamical theory of quantum dot systems

TL;DR

This paper develops a non-equilibrium dynamical theory for quantum dot systems using stochastic fermionic quantum state diffusion, revealing complex Coulomb blockade phenomena and non-Markovian effects in electron transport.

Contribution

It introduces a novel non-equilibrium quantum dynamical framework for quantum dots coupled to reservoirs, capturing transient and steady-state behaviors beyond traditional methods.

Findings

Coulomb blockade observed in spin-degenerate and non-degenerate quantum dots.

Non-monotonic current behavior with Coulomb energy in spin-degenerate case.

Significant non-Markovian effects influencing steady-state currents.

Abstract

We theoretically investigate the non-equilibrium quantum dynamical theory of a quantum dot system coupled to fermionic reservoirs using the recently developed stochastic fermionic quantum state diffusion (FQSD) equation. The exact or approximate dynamical equations associated with the FQSD equation can describe the non-equilibrium quantum transport processes beyond the long-time limit leading to a steady state. We study in details the electron transport of a quantum-dot system coupled to two fermionic environments with different chemical potentials. We report the onset of Coulomb blockade in quantum dots in two distinctive cases: one involving a spin degeneracy one-quantum dot model, and the other a specific spin non-degeneracy two-quantum dot model. While the spin degeneracy case shows that the current in the quantum dot may be blockaded non-monotonically with respect to Coulomb energy, the non-degeneracy case exhibits significant non-Markovian effects, and it enables us to study the relations between initial conditions of the dots and the steady state currents.

Equations

$${\hat{H}}_{\mathrm{tot}}={\hat{H}}_{\mathrm{S}}+{\hat{H}}_{\mathrm{I}}+{\hat{H}}_{\mathrm{E}},$$

$$(1-{\bar{n}}_{jk})|0_{b_{jk}^{\prime }}⟩⟨0_{b_{jk}^{\prime }}|+{\bar{n}}_{jk}|1_{b_{jk}^{\prime }}⟩⟨1_{b_{jk}^{\prime }}|$$

$$\to$$

$${\hat{b}}_{jk}$$

$${\hat{c}}_{jk}$$

$${\hat{b}}_{jk}^{\prime }$$

$${\hat{c}}_{jk}^{\prime }$$

$$\sqrt{1-{\bar{n}}_{jk}}|0_{b_{jk}^{\prime }},0_{c_{jk}^{\prime }}⟩+\sqrt{{\bar{n}}_{jk}}|1_{b_{jk}^{\prime }},1_{c_{jk}^{\prime }}⟩\to |0_{b_{jk}},0_{c_{jk}}⟩.$$

$${\hat{H}}_{\mathrm{tot}}^{\prime }$$

$$\Vert \mathbf{\xi }⟩\equiv \prod _{j=1,2;k}(1-{\xi }_{b_{jk}}{\hat{b}}_{jk}^{†})(1-{\xi }_{c_{jk}}{\hat{c}}_{jk}^{†})|\mathrm{vac}{⟩}_{\mathrm{E}},$$

$$|{\psi }_t({\mathbf{\xi }}^{\ast })⟩$$

$${\hat{\rho }}_r=\int {\mathcal{D}}_g[\mathbf{\xi }]|{\psi }_t({\mathbf{\xi }}^{\ast })⟩⟨{\psi }_t(-\mathbf{\xi })|=\mathcal{M}(|{\psi }_t({\mathbf{\xi }}^{\ast })⟩⟨{\psi }_t(-\mathbf{\xi })|),$$

$${\mathcal{D}}_g[\mathbf{\xi }]\equiv \prod _{\lambda =b,c}\prod _{j=1,2;k}(d{\xi }_{{\lambda }_{jk}}^{\ast }d{\xi }_{{\lambda }_{jk}}{e}^{-{\xi }_{{\lambda }_{jk}}^{\ast }{\xi }_{{\lambda }_{jk}}}),$$

$${\partial }_t|{\psi }_t({\mathbf{\xi }}^{\ast })⟩$$

$${\xi }_{b_{j}t}^{\ast }\equiv -i\sum _k{t}_{b_{jk}}{e}^{i{\omega }_{k}t}{\xi }_{b_{jk}}^{\ast },K_{b_j}(t,s)=\sum _k{t}_{b_{jk}}^2{e}^{-i{\omega }_k(t-s)},$$

$${\xi }_{c_{j}t}^{\ast }\equiv -i\sum _k{t}_{c_{jk}}{e}^{-i{\omega }_{k}t}{\xi }_{c_{jk}}^{\ast },K_{c_j}(t,s)=\sum _k{t}_{c_{jk}}^2{e}^{i{\omega }_k(t-s)},$$

$${\bar{Q}}_{{\lambda }_j}(t,{\mathbf{\xi }}^{\ast })\equiv {\int }_0^{t}ds{K}_{{\lambda }_j}(t,s){\hat{Q}}_{{\lambda }_j}(t,s,{\mathbf{\xi }}^{\ast }),(\lambda =b,c),$$

$${\hat{Q}}_{{\lambda }_j}(t,s,{\mathbf{\xi }}^{\ast })|{\psi }_t({\mathbf{\xi }}^{\ast })⟩\equiv {\overrightarrow{\delta }}_{{\xi }_{{\lambda }_{j}s}^{\ast }}|{\psi }_t({\mathbf{\xi }}^{\ast })⟩,$$

$${\partial }_t{\hat{\rho }}_r$$

$$\displaystyle{\mathcal{M}}\big{(}\bar{Q}_{b_{j}}(t,{\bm{\xi}}^{*})\hat{P}-(\bar{Q}_{c_{j}}(t,{\bm{\xi}}^{*})\hat{P})^{\dagger}\big{)}$$

$$\displaystyle\partial_{t}\hat{\rho}_{r}\!=\!-i[\hat{H}_{\rm S},\hat{\rho}_{r}]\!+\!\sum_{j}\big{(}[F_{1j}\hat{d}\hat{\rho}_{r}\!+\!F_{2j}\hat{\rho}_{r}\hat{d},\hat{d}^{\dagger}]+{\rm H.c.}\big{)}.\quad\quad$$

$${\partial }_t{\hat{\rho }}_r=-i[{\hat{H}}_{\mathrm{S}},{\hat{\rho }}_r]+({\hat{v}}_{1d}-{\hat{v}}_{d2}),$$

$${\hat{v}}_{1d}$$

$${\hat{v}}_{d2}$$

$${\partial }_t\mathrm{Tr}({\hat{\rho }}_r(t)\hat{N})$$

$${\hat{v}}_{1d}$$

$${\hat{v}}_{d2}$$

$$I_{1d}$$

$$I_{d2}$$

$$t_{jk}^2={\Lambda }^2\Omega \Delta \omega \frac{b^2}{({\omega }_k/\Omega -1{)}^2+b^2},$$

$${\hat{H}}_{\mathrm{tot}}$$

$${\xi }_{b_{j\sigma }t}^{\ast }$$

$${\xi }_{c_{j\sigma }t}^{\ast }$$

$${\bar{Q}}_{\nu }(t,{\mathbf{\xi }}^{\ast })\equiv {\int }_0^{t}ds{K}_{\nu }(t,s){\hat{Q}}_{\nu }(t,s,{\mathbf{\xi }}^{\ast }),$$

$${\hat{Q}}_{\nu }(t,s,{\mathbf{\xi }}^{\ast })|{\psi }_t⟩\equiv {\overrightarrow{\delta }}_{{\xi }_{\nu s}^{\ast }}|{\psi }_t({\mathbf{\xi }}^{\ast })⟩,$$

$$K_{b_j}(t,s)$$

Full text

I Introduction

As a potential candidate of quantum devices in the implementation of quantum information processing Loss ; Nori1 ; Nori2 ; Ladd ; Nori3 , quantum dot systems and their quantum transport phenomena have been studied extensively from different points of view Meir ; Wiseman ; Ventra . From a viewpoint of quantum open system, the source and drain in the electron transport may be modelled by two thermal equilibrium reservoirs with two different chemical potentials. Traditionally, treatments of the interaction of a small quantum system with an environment are based on a variety of approaches, such as the many-body Schrödinger equation, the non-equilibrium Green’s function method xx ; Hu0 , the input-output approach InputOutput , the path integral approach hu ; Feynman-Vernon ; Leggett ; Hu1 , and the stochastic Schrödinger equation (SSE) method Gisin ; QuanJump1 ; QuanJump2 ; QuanJump3 ; Diosi1 ; Diosi2 ; Strunz ; Jing . More specifically, the SSE method represented by a non-Markovian quantum state diffusion (QSD) equation can simulate the non-Markovian features of the open quantum systems arising in many interesting scenarios such as a structured environment, a time-delayed external control, or the strong system-environment coupling etc. Such a non-Markovian approach will be useful for studying the temporal behaviour of quantum transport processes in different time scales that is of interest from recent attempts in understanding quantum decoherence and non-equilibrium dynamics. Several approaches have been studies recently including a non-Markovian master equation (ME) for quantum dots coupled to fermionic non-Markovian environments derived by using the Feynman-Vernon influence functional (IF) approach ZhangDQD . A different treatment based on the fermionic non-Markovian stochastic Schrödinger equation (NMSSE) method has been developed Shi2013c ; Zhao2012d ; Chen2013 , and it has been demonstrated that the exact ME can be derived for a quadratic Hamiltonian. For a system with a more generic non-quadratic Hamiltonian, a systematic perturbation is available for numerical calculations.

The purpose of this paper is to study quantum dynamical processes of a class of quantum dot systems coupled to one or more fermionic reservoirs by using the perturbative NMSSE method. We use the NMSSE method Shi2013c ; Zhao2012d ; Chen2013 to derive non-Markovian MEs for the density operators of quantum dot systems coupled to their fermionic environments. We show, in particular, that the NMSSE with quartic interaction terms (Coulomb interaction terms) plays an important role in describing a non-equilibrium process such as the temporal behaviours of the quantum transport processes. We show that the transport processes with the Coulomb interaction can be studied with an approximate ME, and in case of the absence of the Coulomb interaction, an exact ME can be obtained from NMSSE in a straightforward manner. It should be noted that the Coulomb blockade effects are studied for both spin degeneracy and non-degeneracy cases. In the latter case, we report the observation of a correlated noise arising from two separate fermionic reservoirs.

The paper is organized as follows. In Sec. II, we begin with a brief introduction of the NMSSE and a generic non-Markovian ME derived from the NMSSE for the density operator of a single quantum dot coupled to the source and drain modelled by two fermionic reservoirs. Then we discuss the effect of Coulomb blockade in two model systems in Sec. III. In the spin degeneracy one-quantum dot model, we investigate the influence of the bandwidth and chemical potential on the Coulomb blockade effect. In addition, in the spin non-degeneracy two-dot model, with the two-quantum dot resonance condition we show that the steady state current is sensibly dependent on the initial states of the quantum dots. Our discussions and final comments are presented in Sec. IV. Appendix A gives the details of the numerical algorithm for the evaluation of the coefficients of the exact ME in Sec. II.

II Transport Between Two environments through a quantum dot

To begin with, we consider a simple system composed of a single quantum dot coupled to two fermionic environments which correspond to a source and a drain, respectively, as shown in Fig. 1. The environments may represent electrodes or impurities.

We assume that the potential barriers between the quantum dot and environments can localize the electron’s mode in the dot, but still permits the electron tunnel processes. In this case, the Hamiltonian of the total system can be written in the framework of system-plus-environment (we set $\hbar=1$ throughout the paper):

[TABLE]

where $\hat{H}_{\rm S}=\Omega\hat{d}^{\dagger}\hat{d}$ is the Hamiltonian of the quantum dot, $\hat{H}_{\rm E}=\sum_{k}\omega_{k}(\hat{b}^{\prime}_{1k}{}^{\dagger}\hat{b}^{\prime}_{1k}+\hat{b}^{\prime}_{2k}{}^{\dagger}\hat{b}^{\prime}_{2k})$ is the Hamiltonian of the source and drain, and $\hat{H}_{\rm I}=\sum_{k}\sum_{j=1,2}t_{jk}\left(\hat{d}^{\dagger}\hat{b}_{jk}^{\prime}+\hat{b}_{jk}^{\prime}{}^{\dagger}\hat{d}\right)$ is the interaction Hamiltonian between the dot and the environments (source and drain). Here, $\hat{d}^{\dagger}$ and $\hat{d}$ are respectively the fermionic creation and annihilation operators of the electrons in the quantum dot, and $\hat{b}^{\prime}{}^{\dagger}_{jk}$ and $\hat{b}^{\prime}_{jk}$ are respectively the fermionic creation and annihilation operators of the $k$th mode of the $j$th environment ($j=1,2$), and $t_{jk}$ stands for the coupling strength between the quantum dot and the $k$th mode of the $j$th environment. The fermionic operators satisfy the well-known anticommutation relations: $\{\hat{d}^{\dagger},\hat{d}\}=1$, $\{\hat{b}^{\prime}{}^{\dagger}_{jk},\hat{b}^{\prime}_{j^{\prime}k^{\prime}}\}=\delta_{kk^{\prime}}\delta_{jj^{\prime}}$, and $\{\hat{b}^{\prime}{}^{\dagger}_{jk},\hat{d}\}=\{\hat{d}^{\dagger},\hat{b}^{\prime}_{jk}\}=0$. For simplicity, we have assumed that $t_{jk}$ are real numbers without losing generality.

The initial states of the environments are set as two thermal equilibrium states at the same temperature $T$, but with different chemical potentials $\mu_{j}$. The NMSSE method for the thermal environments can be treated by using a Bogoliubov transformation to convert a thermal environment problem to that for a vacuum environment. For this purpose, we first decompose the thermal state into an entangled pure state:

[TABLE]

$\bar{n}_{jk}$ is the average number of electrons in the $k$th mode of the $j$th environment $\omega_{k}\hat{b}^{\prime}_{jk}{}^{\dagger}\hat{b}^{\prime}_{jk}$, which satisfies the Fermi-Dirac distribution $\bar{n}_{jk}=1/(1+e^{\beta(\omega_{k}-\mu_{j})})$, here $\beta$ is the inverse temperature $\beta=1/(k_{B}T)$. $c^{\prime}_{jk}$ is a fictitious mode corresponding to the negative energy: $-\omega_{k}\hat{c}^{\prime}_{jk}{}^{\dagger}\hat{c}^{\prime}_{jk}$ which may be understood as a electron hole. The introduction of those electron-hole pairs can essentially transform the finite-temperature problem into a zero-temperature one. The hole system introduced this way is not interacting with the electron reservoirs, so it does not affect the system dynamics, but simply serves as ancilla qubits for purification of the environments. In fact, if we trace out all the fictitious modes of the entangled pure state, the thermal state of $b^{\prime}_{jk}$ can be obtained. Then we apply a Bogoliubov transformation to convert such an entangled pure state into a vacuum. More specifically, if we choose the following transformation:

[TABLE]

while its inverse is

[TABLE]

then the entangled state in Eq. (2) may be converted to a vacuum corresponding to the annihilation operators $\hat{b}_{jk}$ and $\hat{c}_{jk}$:

[TABLE]

Following the transformation, the Hamiltonian (1) the takes the following form:

[TABLE]

where $t_{b_{jk}}\equiv\sqrt{1-\bar{n}_{jk}}t_{jk}$ and $t_{c_{jk}}\equiv\sqrt{\bar{n}_{jk}}t_{jk}$. Therefore, if we choose the initial state of the environment modes $b_{jk}$ and $c_{jk}$ in the Hamiltonian (6) to be a vacuum state, then by tracing over all $c^{\prime}_{jk}$ modes we will get the Hamiltonian (1) involving only the environment modes $b^{\prime}_{jk}$, which are initially in a thermal state.

In what follows, we will develop an NMSSE approach and define the fermionic quantum trajectory for the quantum dot system. The NMSSE approach will involve the Grassmann-Bargmann (GB) coherent state representation  Glauber1999 , which may stimulate an anti-commutating Gaussian noise. The GB coherent state representation is defined as

[TABLE]

where $\xi_{b_{jk}}$ and $\xi_{c_{jk}}$ are independent Grassmann variables, and ${\bm{\xi}}$ is a collective notation. A fermionic quantum trajectory can be defined as the inner product of the GB coherent states and the wave function of the total system $|\psi_{\rm tot}^{I}(t)\rangle$:

[TABLE]

The reduced density operator of the quantum dot system can be obtained by taking the mean value of the fermionic quantum trajectories over the Grassmann variables:

[TABLE]

where the Grassmann-Gaussian measure is defined as

[TABLE]

and ${\mathcal{M}}(\cdots)$ stands for the fermionic stochastic mean value. The time evolution of the trajectories is governed by the NMSSE:

[TABLE]

where the stochastic processes and Q-operators are defined as follows,

[TABLE]

in which $\overrightarrow{\delta}_{\xi_{\lambda_{j}s}^{*}}$ means the left-functional-derivative with respect to noise $\xi_{\lambda_{j}s}^{*}$, and $K_{\lambda_{j}}(t,s)$ are correlation functions defined as $K_{\lambda_{j}}(t,s)\equiv{\mathcal{M}}(\xi_{\lambda_{j}t}\xi_{\lambda_{j}s}^{*})$. Eq. (11) is also called the fermionic quantum state diffusion (FQSD) equation.

When we apply an extended Novikov theorem (for the bosonic case, see Ref. YDGS99 ) to calculate the stochastic mean value: {\mathcal{M}}\big{(}\xi^{*}_{\lambda_{j}t}\hat{P}\big{)}, where $\hat{P}\equiv|\psi_{t}({\bm{\xi}}^{*})\rangle\langle\psi_{t}(-{\bm{\xi}})|$, a formal non-Markovian ME is derived:

[TABLE]

Using the Heisenberg approach, the mean value terms above can be exactly evaluated:

[TABLE]

where the expressions of $F_{ij}$ can be found in Appendix A. By substituting Eq. (14) into Eq. (13), we obtain the ME:

[TABLE]

Analogue to the classical continuity equation: $\partial_{t}\rho+\nabla\cdot\vec{j}=0$, the current of the quantum dot can be connected to the change rate of the particle number. We split the formal ME (13) into a system part and two environment-involved parts:

[TABLE]

where $\hat{v}_{1d}$ and $\hat{v}_{d2}$ are defined as

[TABLE]

Usually $\hat{H}_{\rm S}$ commutes with the system particle number operator $\hat{N}$, so we have

[TABLE]

Corresponding to the $\nabla\cdot\vec{j}$ term, we can define I_{1d}\equiv q_{e}{\rm Tr}\big{(}\hat{v}_{1d}\hat{N}\big{)}, and I_{d2}\equiv q_{e}{\rm Tr}\big{(}\hat{v}_{d2}\hat{N}\big{)}, which stand for the expectation values of currents with respect to directions: $1\rightarrow d$ and $d\rightarrow 2$ respectively. Here $1$ and $2$ denote the $1$-st and the $2$-nd environment, $d$ denotes the quantum dot, and $q_{e}$ is the charge of one electron.

By comparing Eq. (14) with Eq. (17), the currents may be obtained:

[TABLE]

where the elements $\rho_{00}$ and $\rho_{11}$ are defined as $\hat{\rho}_{r}\equiv\rho_{00}|0\rangle\langle 0|+\rho_{11}|1\rangle\langle 1|$, and the superscript $R$ denotes the real part of a complex number.

For simplicity, we use the Lorentzian spectrum in our numerical simulations:

[TABLE]

where $\varDelta\omega$ is the eigen-frequency spacing of the discrete spectrum, the two environments are represented by the same relative bandwidth $b$ and the same coupling strength factor $\Lambda$. In Fig. 2, we can see how the bandwidth $b$ influence the currents in the extreme non-Markovian case. When $b$ decreases, the amount of the currents reduces to [math], and the oscillation start to take place. The difference between $I_{1d}$ and $I_{d2}$ is related to {\rm Tr}\big{(}\hat{d}^{\dagger}\hat{d}\,\partial_{t}\hat{\rho}_{r}(t)\big{)}, so if $\hat{\rho}_{r}$ has a steady state, $\varDelta I\equiv I_{1d}-I_{d2}$ will converge to [math]. In Fig. 3, we can see that the average current: $I\equiv(I_{1d}+I_{d2})/2$ is smoother than $I_{1d}$ and $I_{d2}$,

and display a simpler dissipated oscillation pattern. Actually, in this model $I(t)$ will not be affected by $\hat{\rho}_{r}(0)$ due to the symmetry between the spectrum $t_{1k}$ and $t_{2k}$.

In Fig. 3, we can see that when the bandwidth

becomes bigger, the current $I$ approaches the steady state faster, the oscillation of the current gets weaker, and beyond a threshold of the bandwidth the oscillation disappears. To evaluate how fast the current $I$ approaches its steady state in the extreme non-Markovian case, we plot the stationary point of the current (where the time derivative equals 0) in Fig. 4.

Fig. 4 shows that the chemical potential and temperature have no visible influence on the value $t_{p}$, where $t_{p}$ is the time coordinate corresponding to the first turning point of the current. Therefore, it shows that the bandwidth $b$ is the only important factor that can help $I$ to reach its steady value faster. In Fig. 5, the irrelevances between $\mu_{1}$, $T$ and $t_{p}$ are further demonstrated and we can see how the bandwidth influences the $t_{p}$. In both (a) and (b), the bandwidth range starts from $0.001$.

Now we use ME to discuss the steady state current denoted as $I_{st}\equiv I(t\!\rightarrow\!\infty)$. In Fig. 6, we plot the 3-D figure of $I_{st}(b,\mu_{1})$, and we see that there are two relative flat regions corresponding to the top part and bottom part respectively, and between these two there are two relative linear regions. Note that the bandwidth $b$ starts from $0.002$ in Fig. 6.

III Effect of the Coulomb blockade in non-Markovian environments

We now turn to the problem of the transport involving the Coulomb blockade effect. When there are more than one electrons in the dot or there are more than one dots in the system, the Coulomb repulsion between electrons or dots may blockade the current Gurvitz ; you1999 ; Yan1 ; Yan2 . Such an effect can be applied to measure the state of the dot or control the current. In this section, we will use approximate MEs to investigate the Coulomb blockade effect in the non-Markovian case. We consider two separate cases. In the first case, the spin degrees of freedom will be taken into account, the second case will address a non-degeneracy situation, which will lead to interesting correlated noises phenomena.

III.1 Coulomb blockade in the spin degeneracy situation

In this model we suppose all the energy levels in the environments or in the dot contain two modes which correspond to spin up and spin down, respectively. We assume that when the electrons tunnel in or tunnel out of the dot, their spins do not flip, as illustrated in Fig. 7. In this case, the Hamiltonian can be written as

[TABLE]

where $\hat{H}_{\rm S}=\sum_{\sigma=\downarrow,\uparrow}\Omega\hat{d}_{\sigma}^{\dagger}\hat{d}_{\sigma}+\Omega_{c}\hat{d}_{\downarrow}^{\dagger}\hat{d}_{\downarrow}\hat{d}_{\uparrow}^{\dagger}\hat{d}_{\uparrow}$ is the Hamiltonian of the quantum dot, and $\Omega_{c}$ is the the Coulomb repulsion energy between the electrons.

The method used here closely follows the presentation in Sec. II. By applying the purification and Bogoliubov transformations to the two environments in thermal equilibrium states, we introduce eight independent Grassmann Gaussian stochastic processes:

[TABLE]

and altogether sixteen Q-operators corresponding to these Grassmann stochastic processes:

[TABLE]

where the non-commutative noises satisfy ${\mathcal{M}}(\xi^{*}_{\nu t})={\mathcal{M}}(\xi^{*}_{\nu t}\xi^{*}_{\mu^{\prime}s})=0$, and ${\mathcal{M}}(\xi_{\nu t}\xi^{*}_{\nu^{\prime}s})=\delta_{\nu\nu^{\prime}}K_{\nu}(t,s)$. Here $K_{\nu}(t,s)$ are the correlation functions of the noises, and the index $\nu=b_{1\downarrow},b_{1\uparrow},b_{2\downarrow},b_{2\uparrow},c_{1\downarrow},c_{1\uparrow},c_{2\downarrow},c_{2\uparrow}$. Since the spectrum is symmetric with respect to spin interchange, the correlation functions satisfy $K_{\lambda_{j\downarrow}}(t,s)=K_{\lambda_{j\uparrow}}(t,s)$, where $\lambda=b,c$ and $j=1,2$. Thus, we will eliminate the spin index in the correlation functions in this subsection hereinafter:

[TABLE]

Then the fermionic QSD equation can be derived:

[TABLE]

We can see that the functional form of the coupling terms with $\xi^{*}_{\lambda_{1\sigma}t}$ in Eq. (25) is the same as the ones with $\xi^{*}_{\lambda_{2\sigma}t}$, so their functional derivatives should equal: $\overrightarrow{\delta}_{\xi^{*}_{\lambda_{1\sigma}s}}|\psi_{t}({\bm{\xi}}^{*})\rangle=\overrightarrow{\delta}_{\xi^{*}_{\lambda_{2\sigma}s}}|\psi_{t}({\bm{\xi}}^{*})\rangle$, which means $\hat{Q}_{\lambda_{1\sigma}}=\hat{Q}_{\lambda_{2\sigma}}$. Thus, we will eliminate the index $j$ in $\hat{Q}_{\lambda_{j\sigma}}$ in this subsection hereinafter. $\bar{Q}_{\nu}$ can also be simplified. We define $K_{\lambda}(t,s)\equiv\sum_{j}K_{\lambda_{j}}(t,s)$, and $\bar{Q}_{\lambda_{\sigma}}\equiv\sum_{j}\bar{Q}_{\lambda_{j\sigma}}$, then we have $\bar{Q}_{\lambda_{\sigma}}(t,{\bm{\xi}}^{*})=\int_{0}^{t}dsK_{\lambda}(t,s)\hat{Q}_{\lambda_{\sigma}}(t,s,{\bm{\xi}}^{*})$.

The evolution of the Q-operators is governed by the equations and boundary conditions:

[TABLE]

which can be derived from the consistency conditions: $\partial_{t}\overrightarrow{\delta}_{\!\xi^{*}_{\lambda_{j\sigma}s}}|\psi_{t}({\bm{\xi}}^{*})\rangle=\overrightarrow{\delta}_{\!\xi^{*}_{\lambda_{j\sigma}s}}\partial_{t}|\psi_{t}({\bm{\xi}}^{*})\rangle$. In this model, the Q-operators will involve the noise terms up to infinite degree, however when the coupling strength between the system and environments is weak, the zeroth order approximation is proved to be a good approximation for perturbation.

If we neglect all the noise terms in Eq. (27), an approximate equation can be obtained:

[TABLE]

We use $\hat{Q}_{\lambda_{\sigma}}^{\{0\}}$ and $\bar{Q}_{\lambda_{\sigma}}^{\{0\}}$ to denote these approximate Q-operators. Then Eq. (28) can be rewritten as

[TABLE]

where the boundary conditions are $\hat{Q}_{b_{\sigma}}^{\{0\}}(t=s,s)=\hat{d}_{\sigma}$ and $\hat{Q}_{c_{\sigma}}^{\{0\}}(t=s,s)=-\hat{d}_{\sigma}^{\dagger}$. When the Coulomb interaction is present, the zeroth order Q-operators are of such forms:

[TABLE]

Due to the symmetry with respect to spin interchange in Eq. (29), we have the following relations: $f_{\lambda_{\downarrow}}=f_{\lambda_{\uparrow}}$, $g_{\lambda_{\downarrow}}=g_{\lambda_{\uparrow}}$, $F_{\lambda_{j\downarrow}}=F_{\lambda_{j\uparrow}}$, $G_{\lambda_{j\downarrow}}=G_{\lambda_{j\uparrow}}$, $\lambda=b,c$, and $j=1,2$. Subsequently, we may drop the index $\sigma$ in $f_{\lambda_{\sigma}}$, $g_{\lambda_{\sigma}}$, $F_{\lambda_{j\sigma}}$ and $G_{\lambda_{j\sigma}}$ in this subsection from now on. By substituting the expansions in Eq. (30) into Eq. (29), the equations of $f_{\lambda}$ and $g_{\lambda}$ can be obtained:

[TABLE]

where the boundary conditions are $f_{b}(t\!=\!s,s)\!=\!1$, $f_{e}(t\!=\!s,s)\!=\!-1$, $g_{\lambda}(t\!=\!s,s)\!=\!0$.

Once the coefficients of Q-operators are solved, we then can easily derive the ME for the density operator of the system:

[TABLE]

where $F_{\lambda}\equiv\sum_{j=1,2}F_{\lambda_{j}}$, and $G_{\lambda}\equiv\sum_{j=1,2}G_{\lambda_{j}}$; the superscript $R$ stands for the real part of a complex number.

Following the same procedure introduced in Sec. II, we can find the electronic current between the quantum dot and the $1$st environment:

[TABLE]

where we expand the reduced density operator as

[TABLE]

The current $I_{d2}$ can be obtained similarly:

[TABLE]

In the numerical simulation part, we still use the Lorentzian spectrum:

[TABLE]

where for simplicity the couplings to the source and drain are assumed to be symmetric.

To make the approximation valid, we set $\Lambda=0.014$ and $b\geq 0.05$ in this subsection. In this model the steady state of $\hat{\rho}_{r}$ can be fully determined by the coefficients $F_{\lambda}$ and $G_{\lambda}$, thus we do not indicate the initial state of the system. From Fig. 8, we can see how the Coulomb energy $\Omega_{c}$ influences the steady state current $I_{st}\equiv(I_{1d}(t\!\rightarrow\!\infty)+I_{d2}(t\!\rightarrow\!\infty))/2$. In all the four subfigures which represent different bandwidths, there is a clear line corresponding to $\mu_{1}=\Omega+\Omega_{c}$. In the region of $\mu_{1}\leq\Omega+\Omega_{c}$, the steady state current keeps almost a constant value which is just the blockaded current. Fig. 9 displays some sections of Fig. 8 perpendicular to axis $\mu_{1}$. We can see that, when the Coulomb interaction is ignored, the steady state current is given by $-$9 nA, when $\Omega_{c}\rightarrow\infty$, the steady state current converges to about $2/3$ of $-$9 nA. Before reaching $2/3$ of the current $I_{st}(\Omega_{c}=0)$, the steady state current may have higher values, which are dependent on $\mu_{1}$. When $\mu_{1}$ increases, the overshoot will converge to $1/2$ of the current $I_{st}(\Omega_{c}=0)$.

III.2 Coulomb blockade in the correlated noise situation

Suppose two quantum dots are arranged in parallel between the source and drain, and the spin degrees of freedom of all electrons are neglected, (the non-degeneracy case can be achieved by adding a strong magnetic field) then for each electron in the environments, it has a chance to tunnel into either dot. A schematic configuration for this situation is shown in Fig. 10, as to be shown below, this case will generate correlated noises.

The Hamiltonian of the total system may be written as

[TABLE]

where $\hat{H}_{\rm S}=\sum_{\epsilon=1,2}\Omega_{\epsilon}\hat{d}_{\epsilon}^{\dagger}\hat{d}_{\epsilon}+\Omega_{c}\hat{d}_{1}^{\dagger}\hat{d}_{1}\hat{d}_{2}^{\dagger}\hat{d}_{2}$ is the Hamiltonian of the dots, and $\Omega_{c}$ is the Coulomb repulsive energy. Again, the parameter $t_{\epsilon,jk}$ is taken as a real number representing the coupling between the $\epsilon$th dot to the $k$th mode of the $j$th environment ($j=1,2$).

A set of Grassmann stochastic processes can be obtained by applying the purification and Bogoliubov transformations to the environments,

[TABLE]

which satisfy ${\mathcal{M}}(\xi^{*}_{\epsilon\lambda_{j}t})={\mathcal{M}}(\xi^{*}_{\epsilon\lambda_{j}t}\xi^{*}_{\epsilon^{\prime}\lambda^{\prime}_{j^{\prime}}s})=0$. The noise $\xi^{*}_{\epsilon\lambda_{j}t}$ represents the stochastic influence of the environment $\lambda_{j}$ on the dot $\epsilon$. Usually $\xi^{*}_{\epsilon\lambda_{j}t}$ and $\xi^{*}_{\epsilon^{\prime}\lambda_{j}s}$ are highly correlated, the elements of the cross correlation function matrix is defined as

[TABLE]

where the subscript $\nu$ is a short notation of $\lambda_{j}$ which represents the four possible indices $b_{1}$, $b_{2}$, $c_{1}$, $c_{2}$.

Then the equation of motion for the fermionic quantum trajectory can be obtained:

[TABLE]

The functional derivatives in Eq. (40) can be replaced by Q-operators:

[TABLE]

Then the fermionic QSD equation is finally obtained:

[TABLE]

By a similar observation on the coupling forms of $\xi^{*}_{\epsilon\lambda_{1}t}$ and $\xi^{*}_{\epsilon\lambda_{2}t}$ in Eq. (42), we conclude that the Q-operators should satisfy: $\hat{Q}_{\epsilon\lambda_{1}}=\hat{Q}_{\epsilon\lambda_{2}}$. Thus, we may drop the index $j$ in $\hat{Q}_{\epsilon\lambda_{j}}$ from now on in this subsection.

When the couping is weak, we may use $\hat{Q}_{\epsilon\lambda}^{\{0\}}$ to replace the exact Q-operators, where the approximate operators evolve according yo the following equations:

[TABLE]

We expand $\bar{Q}_{\epsilon\lambda_{j}}^{\{0\}}$ and $\hat{Q}_{\epsilon\lambda}^{\{0\}}$ with respect to the annihilation and creation operators of the dots:

[TABLE]

By substituting these expansions into Eq. (43), we get the equations and boundary conditions of the coefficients:

[TABLE]

where the explicit expressions of $A$, $B$, $C$, $D$ and $E$ are

[TABLE]

When the coefficients above are solved, we directly obtain the zeroth order ME through the extended Novikov theorem:

[TABLE]

Then the current $I_{1d}$ and $I_{d2}$ can be separated from the ME:

[TABLE]

where the reduced density operator is expanded as

[TABLE]

When the resonance condition is satisfied $\Omega_{1}=\Omega_{2}$, it is easy to see that the steady state current will be affected by the initial condition of $\hat{\rho}_{r}$, which is different from the model in the previous subsection. In numerical simulations, we use such a Lorentzian spectrum form:

[TABLE]

With the chosen spectrum, when the resonance condition is valid ($\Omega_{1}=\Omega_{2}$), then all the four components in the cross correlation function are identical: $C_{11,\lambda_{j}}=C_{12,\lambda_{j}}=C_{21,\lambda_{j}}=C_{22,\lambda_{j}}$. By substituting these relations into Eq. (45) and (46), we get $f_{11\lambda}=f_{12\lambda}=f_{21\lambda}=f_{22\lambda}$, $g_{11\lambda}=g_{12\lambda}=g_{21\lambda}=g_{22\lambda}$, $F_{11\lambda_{j}}=F_{12\lambda_{j}}=F_{21\lambda_{j}}=F_{22\lambda_{j}}$ and $G_{11\lambda_{j}}=G_{12\lambda_{j}}=G_{21\lambda_{j}}=G_{22\lambda_{j}}$. To simplify the notation, we drop the first two indices of the above coefficients. With these symmetries, Eq. (LABEL:Coulomb2current1) can be simplified to

[TABLE]

Due to the restriction on fermionic density operators, we only need to consider some of the matrix elements, which can be rearranged into a row vector, then the ME may be written as

[TABLE]

where the superscript $T$ denotes the transpose of a matrix, and $M$-matrix is of the form:

[TABLE]

The explicit expressions of the coefficients in $M$ are

[TABLE]

where the superscript $I$ stands for the imaginary part of a complex number.

The rank of $M$ contains a wealth of information about the steady-state Coulomb blockade. $M$ is generally a rank four matrix, however, the rank can be five when the resonance condition is satisfied, which reduces to the degeneracy case considered in the previous subsection. Therefore, we will be focused on the rank four case where the steady state current $I_{st}\!\equiv\!(I_{1d}(t\!\rightarrow\!\infty)+I_{d2}(t\!\rightarrow\!\infty))/2$ is shown to be sensitively dependent on the initial state $\hat{\rho}_{r}(0)$. Clearly, by the linearity of MEs, we only to consider the steady states of a complete basis of $\hat{\rho}_{r}(0)$, which can generate an arbitrary initial $\hat{\rho}_{r}(0)$ by a convex combination. For a reason to become clear later, we select following eight specific initial states as the basis ($i=1,...,8$):

[TABLE]

Obviously, the choice of a basis is not unique. In Fig. 11, we show that the eight initial basis states may be divided into three groups. The group I consists of first two initial vectors in Eq. (55), Group II contains four initial states ranging from 3 to 6; Group III contains the rest two initial states (7,8). As seen from Fig. 11, the initial states in each group will lead to the identical steady states, that is, $I_{st,1}=I_{st,2}$, $I_{st,3}=I_{st,4}=I_{st,5}=I_{st,6}$, and $I_{st,7}=I_{st,8}$, where the index $i$ in $I_{st,i}$ denotes the $i$th initial state of $\hat{\rho}_{r}$ in Eq. (55).

From Fig. 11, we can see that when $\Omega_{c}=0$ all the initial states lead to the same steady state current. When $\Omega_{c}\neq 0$, there is no Coulomb blockade effect in Group I. However, it is interesting to notice that in Group III, $I_{st}$ is blockaded monotonically with respect to $\Omega_{c}$, and it can be blockaded to [math] when $\Omega_{c}$ is high enough. In Group II, $I_{st}$ is also blockaded monotonically, however it can only be blockaded to near half of the original value $I_{st}(\Omega_{c}=0)$.

In order to get a more detailed picture about the Coulomb blockage for different coulomb coupling energies, Figs. 12 and 13 are plotted to show how $\Omega_{c}$ influences $I_{st}$ with different chemical potential $\mu_{1}$ and the bandwidth $b$. Different from the spin degeneracy case (see Fig. 9), the blockade is always monotonic with respect to $\Omega_{c}$. Particularly, we notice the similarity between the curves in Fig. 12 and curves in Fig. 13. In fact, we can show that the currents $I_{st,3}$ and $I_{st,1}+I_{st,7}$ in Fig. 14 are approximately equal, that is, $2I_{st,3}=I_{st,1}+I_{st,7}$.

More theoretical analysis of three groups of initial states can shed light on the steady state classifications. From Eq. (51) we see that due to the symmetry of the resonance the currents $I_{1d}$ and $I_{d2}$ only involve some elements of $\hat{\rho}_{r}$: $\rho_{00}$, $(\rho_{11}+\rho_{22})/2$, $\rho_{33}$, and $\rho_{12}^{R}$. According to Eq. (52) and (53), a set of closed equations can be found for these elements:

[TABLE]

where $\rho_{+}\equiv(\rho_{11}+\rho_{22})/2$. Obviously, from Eq. (51) and (56), we can derive: $I_{1d,3}(t)=I_{1d,4}(t)=I_{1d,5}(t)=I_{1d,6}(t)$, and $I_{d2,3}(t)=I_{d2,4}(t)=I_{d2,5}(t)=I_{d2,6}(t)$, where $I_{1d,i}$ and $I_{d2,i}$ $(i=1,...,8)$ denote the currents corresponding to the eight initial states mentioned in Eq. (55). Eq. (56) can be separated into two parts:

[TABLE]

When $i=7,8$, the initial states satisfy $\rho_{00}(0)=\rho_{+}(0)+\rho_{12}^{R}(0)=0$, then according to the linearity of Eq. (58), we have $\rho_{00}(t)=\rho_{+}(t)+\rho_{12}^{R}(t)=0$. Thus, Eq. (57) turns into

[TABLE]

and the trace formula of the $\hat{\rho}_{r}(t)$ turns into $2\rho_{+}+\rho_{33}=1$. By combining Eq. (59) and the trace formula, we can conclude that the four elements of $\hat{\rho}_{r}(t)$ in Eq. (56) corresponding to the 7-th and 8-th initial states will converge to the same steady states: $\rho_{+}(t\!\rightarrow\!\infty)=-\rho_{12}^{R}(t\!\rightarrow\!\infty)=F_{st}/(2(F_{st}-C_{st}))$, $\rho_{33}(t\!\rightarrow\!\infty)=C_{st}/(C_{st}-F_{st})$, and $\rho_{00}(t)=0$, where the subscript $st$ denotes the steady value of a coefficient. Thus, we have $I_{st,7}=I_{st,8}$.

A similar analysis can be done for Group I with $i=1,2$. Then $\rho_{33}(t)=0$, and $\rho_{+}(t)=\rho_{12}^{R}(t)$ will be obtained. The steady states in these cases are given by: $\rho_{00}(t\rightarrow\infty)=B_{st}/(B_{st}-A_{st})$, $\rho_{+}(t\rightarrow\infty)=\rho_{12}^{R}(t\rightarrow\infty)=A_{st}/(2(A_{st}-B_{st}))$. So we have $I_{st,1}=I_{st,2}$, and $\rho_{33}(t)=0$. Since $\rho_{33}$ remains zero for the entire evolution, hence the Coulomb blockade effect does not take place in Group I.

In Eq. (55), we can see that $\hat{\rho}_{r,3}(0)+\hat{\rho}_{r,4}(0)=\hat{\rho}_{r,2}(0)+\hat{\rho}_{r,8}(0)$. Due to the linearity of the ME and the expressions in Eq. (51), we can conclude that $I_{st,3}+I_{st,4}=I_{st,2}+I_{st,8}$. Then, by combining the conclusions about the three groups of the initial states, the relation $2I_{st,3}=I_{st,1}+I_{st,7}$ can be obtained.

IV conclusions

We have applied the NMSSE approach to derive the ME for the density operator of a quantum dot coupled to its fermionic environments. We have also derived an approximate ME valid to the zeroth order perturbation. We emphasize that the zeroth order perturbation effectively ignores the noise terms, but still retains higher order coupling strength terms, hence it render the perturbation method a very useful tool in describing a weakly non-Markovian dynamics. A higher order noise perturbation beyond the zeroth order approximation is possible if a highly non-Markovian environment is involved and more computational resources are available.

With both the exact and approximate methods, the current between the quantum dots and the source and drain has been discussed. The effect of the Coulomb interaction between the electrons has been discussed on the transport between the quantum dots and the source and drain. We have found that in the spin degeneracy case when $\mu_{1}$ is sufficiently large, the steady state current may not be monotonically blockaded with respect to the Coulomb energy. In a specific correlated noise case, we have found that the steady state current may sensitively depend on the initial state of the quantum dots. More specifically, we have shown that there is a class of initial states (Group I states) for those states the Coulomb blockade does not occur, and for Group III states the maximal Coulomb blockade is observed. All the other initial states can be described by the combination of Group I and Group II states.

acknowledgements

We acknowledge the grant support from the NSF PHY-0925174, the National Natural Science Foundation of China (Grant No. 11304024). We thanks Mr. Quanzhen Ding for useful discussions on the manuscript.

Appendix A The algorithm for the evaluation of the coefficients $F_{ij}$ of

the exact ME

The coefficients $F_{ij}$ are defined as

[TABLE]

where $u_{i}^{\lambda}$ are the coefficients in the expansions of ${\mathcal{M}}(\hat{Q}_{\lambda}\hat{P})$:

[TABLE]

When the correlation functions satisfy $K_{\lambda_{j}}(t+\tau,s+\tau)=K_{\lambda_{j}}(t,s)$, the coefficients $F_{ij}$ can be constructed via the basic solution $u$. Here we only list the main results without showing too many details. We define the correlation functions: $K_{\lambda_{j}}(t-s)\equiv K_{\lambda_{j}}(t,s)$, $K_{\lambda}(t)\equiv\sum_{j=1,2}K_{\lambda_{j}}(t)$, $K(t)\equiv K_{b}(t)+K^{*}_{c}(t)$. $u$ is the solution of the integro-differential equation:

[TABLE]

Then $F_{ij}$ can be constructed through the algorithm:

[TABLE]

$C_{ij}$ can be generated from

[TABLE]

$A_{j}$ and $B_{j}$ can be obtained from solving equations:

[TABLE]

$D_{1}$ and $D_{2}$ can be obtained from solving equations:

[TABLE]

When the model contains lots of symmetries, the current can be obtained even without solving the $\hat{\rho}_{r}$. According to Eq. (19), the average current can be obtained from

[TABLE]

By utilizing the trace formula of the reduced density operator: $\rho_{00}+\rho_{11}=1$, Eq. (67) can be rewritten as

[TABLE]

By substituting Eq. (LABEL:4Apd2F) into the above equation, we get

[TABLE]

And we know that $K_{b_{j}}+K^{*}_{c_{j}}$ is equal to the correlation function in the vacuum case, explicitly for each $j$, we have

[TABLE]

Thus, when $t_{1k}=t_{2k}$, we have $K_{b_{1}}+K^{*}_{c_{1}}=K_{b_{2}}+K^{*}_{c_{2}}$. In such a case, $F_{12}-F_{11}-F_{22}+F_{21}=0$, and the expression of the average current in Eq. (68) can be simplified:

[TABLE]