Non-Perturbative Electron Counting Statistics

TL;DR

This paper introduces a non-perturbative method using hierarchical equations of motion to accurately analyze electron transport statistics in nanoscale systems, especially where traditional perturbative methods fail.

Contribution

It develops a novel non-perturbative framework for full counting statistics applicable to systems with strong system-bath couplings, improving accuracy over perturbative approaches.

Findings

Perturbative methods are inadequate for certain electron transfer regimes.

The new framework accurately captures current fluctuations in complex systems.

Comparison shows significant differences from traditional approaches.

Abstract

The theory of full counting statistics allows complete characterization of charge transport processes through nanoscale systems. The majority of existing theoretical treatments used to obtain the current cumulants rely on perturbative approximations with respect to either the system-bath coupling or the electronic tunnelling coupling within the system. This is not generally well suited to electron transfer through organic molecules where these couplings can be on the same order of magnitude. Here we present a non-perturbative approach that uses the formalism of the hierarchical equations of motion to describe the system-bath dynamics while retaining a weak coupling with respect to source and drain leads. The scope of this new framework is demonstrated by comparisons with a perturbative approach for a dimer system coupled to a thermal bath. The inadequacy of the perturbative approach to describe current fluctuations in a variety of regimes is shown.

Equations

$$H_S$$

$$H_B$$

$$H_L$$

$$H_{SB}$$

$$H_{SL}$$

$$\frac{d}{dt}\chi(t)=-i[H_{S}+H_{B}+H_{SB},\chi(t)]+\mathcal{D}_{L}\big{(}\chi(t)\big{)}$$

$$\mathcal{D}_{L}\big{(}\chi(t)\big{)}=\sum_{j=L,R}\Gamma_{j}\Big{[}s_{j}\rho_{S}s_{j}^{\dagger}-\frac{1}{2}\{s_{j}^{\dagger}s_{j},\rho_{S}\}\Big{]}.$$

$${\dot{\sigma }}^{(\mathbf{n})}(t)=$$

$$\dot{\sigma }(t)=\mathcal{H}\sigma (t)$$

$$\mathcal{H}(\chi )|0(\chi )⟩⟩={\lambda }_0(\chi )|0(\chi )⟩⟩.$$

$${\lambda }_0(\chi )=⟨⟨\tilde{0}|\Delta \mathcal{H}(\chi )|0⟩⟩$$

$$|0(\chi )⟩⟩=|0⟩⟩+\mathcal{Q}|0(\chi )⟩⟩.$$

$$\mathcal{Q}\mathcal{L}(\chi )\mathcal{Q}|0(\chi )⟩⟩=[{\lambda }_0(\chi )-\Delta \mathcal{H}(\chi )]|0(\chi )⟩⟩.$$

$$\mathcal{Q}|0(\chi )⟩⟩=\mathcal{R}[{\lambda }_0(\chi )-\Delta \mathcal{H}(\chi )]|0(\chi )⟩⟩,$$

$$|0(\chi )⟩⟩=|0⟩⟩+\mathcal{R}[{\lambda }_0(\chi )-\Delta \mathcal{H}(\chi )]|0(\chi )⟩⟩.$$

$${\lambda }_0(\chi )$$

$$|0(\chi )⟩⟩$$

$$\Delta \mathcal{H}(\chi )$$

$$⟨⟨I^1⟩⟩=⟨⟨\tilde{0}|{\mathcal{H}}^{(1)}|0⟩⟩,$$

$$⟨⟨I^2⟩⟩=⟨⟨\tilde{0}|{\mathcal{H}}^{(2)}-2{\mathcal{H}}^{(1)}\mathcal{R}{\mathcal{H}}^{(1)}|0⟩⟩$$

$$\mathcal{L}=\left(\begin{array}{ccccc}-{\Gamma }_L& 0& {\Gamma }_R{e}^{i\chi }& 0& 0\\ {\Gamma }_L& 0& 0& 0& 2T_c\\ 0& 0& -{\Gamma }_R& 0& -2T_c\\ 0& {\gamma }_{+}& -{\gamma }_{-}& -\frac{{\Gamma }_R}{2}-\gamma & -ϵ\\ 0& -T_c& T_c& ϵ& -\frac{{\Gamma }_R}{2}-\gamma \end{array}\right)$$

$$\gamma =$$

$${\gamma }_{\pm }=$$

$$F^{(2)}(0)=1+\frac{2}{\frac{1}{4}{\Gamma }_R^2{\Gamma }_L+2T_c^2{\Gamma }_L+4{ϵ}^2{\Gamma }_L+T_c^2{\Gamma }_R}$$

$$H_S=$$

$$\mathcal{D}(\rho (t))=$$

Full text

††preprint: APS/123-QED

Molecular electron transport is a process fundamental to life, being involved in key stages of both photosynthesis and respiration Gray and Winkler (1996). Incorporating DNA strands Fink and Schönenberger (1999) and individual protein units Gerster et al. (2012) into electric circuits has opened up the possibility of investigating electron transport through organic molecules by analysing the current passing through the system. There is also great potential to use molecular electron transport in novel nanoscale electronic devices Tao (2006); Selzer and Allara (2006). The use of single molecule junctions in diodes, switches and transistors can generate a much richer behaviour than similar semiconductor based devices in part due to the strong interaction between electronic and vibrational components of the molecules Tao (2006). This poses a challenging theoretical problem as the electron-vibrational interaction strength must be treated on the same level as the electronic tunnelling coupling in order to accurately describe electron transport.

The full counting statistics of electron transport provides dynamical information about electron correlations encoded in the probability distribution of the number of charges transferred through the system Levitov and Lesovik (1993); Blanter and Büttiker (2000); Marcos et al. (2010). Features of the electron dynamics such as quantum coherence Kießlich et al. (2007) or the system energy level structure Belzig (2005) can be probed with such a technique, while signatures of the vibrational environment may also be present in the current statistics Haupt et al. (2006). Existing frameworks which provide access to the full counting statistics usually involve a perturbative description of either the electron-vibrational coupling or the electronic tunnelling coupling within the system Brandes (2005). In non-Markovian approaches Braggio et al. (2009a); Flindt et al. (2008, 2010) the current cumulants can be derived from a generalized master equation where the system-environment coupling is often not treated exactly and so the parameter regimes in which they are valid are restricted. Some frameworks do allow for arbitrary system-bath coupling while still considering inexact bath dynamics Braggio et al. (2009b). However, non-perturbative frameworks incorporating an exact bath dynamics have been developed for the full counting statistics of bath observables Cerillo et al. (2016).

In the following we report a non-perturbative approach to electron counting statistics which provides an exact treatment of the coupling between the system and environment as well as the electronic tunnelling coupling within the system. Using the hierarchical equations of motion Tanimura (2006); Ishizaki and Fleming (2009); Tanimura (2012); Kreisbeck et al. (2011) we derive an expression for the zero-frequency current cumulants for non-Markovian quantum dissipative systems which has an analogous formulation in the well studied Markovian regime Flindt et al. (2005); Marcos et al. (2010); Flindt et al. (2010). We demonstrate the wide range of applicability of our approach by comparison with a perturbative framework.

A representation of the model we use is shown in Fig. 1. We consider a dimer system coupled to fermionic leads which either supply or take away electrons from the system. The molecular dimer consists of discrete electronic energy levels coupled to a bosonic phonon bath. The Hamiltonian for the full system is $H=H_{S}+H_{B}+H_{L}+H_{SB}+H_{SL}$ where

[TABLE]

are the Hamiltonians for the system, bath and leads respectively, and we have set $\hbar=1$. Here, $\epsilon$ is the energy bias between the left $|L\rangle$ and right $|R\rangle$ dimer sites and $T_{c}$ is the electronic tunnelling coupling between the sites. $\omega_{q}$ is the frequency of the qth bosonic mode while $b^{\dagger}_{q}$ and $b_{q}$ create and annihilate a phonon in the qth mode respectively. Finally $\epsilon_{j,k}$ is the energy of an electron in the kth energy level of lead $j$ while $c^{\dagger}_{j,k}$ and $c_{j,k}$ create and annihilate an electron in the kth energy level of lead $j$ respectively. The system-bath and system-lead coupling are given by

[TABLE]

The electron-phonon coupling between site $j$ and phonon mode $q$ is given by $g_{j,q}$ while $\nu_{j,k}$ is the coupling between site $j$ and kth energy level in lead $j$. The operators $s^{\dagger}_{j}$ and $s_{j}$ create or annihilate an electron on site $j$ respectively. The coupling between the dimer and the leads is assumed to be weak such that a second order perturbative expansion can be carried out along with the Born-Markov approximation to obtain the action of the leads on the reduced system Brandes (2005); Harbola et al. (2006). A further assumption that the electronic energy levels lie well within the voltage bias window between the leads is made so the infinite bias limit can be applied Harbola et al. (2006). This gives the equation of motion for the system and bath in the Schrödinger picture

[TABLE]

where $\chi(t)=tr_{L}\{\rho_{SB}(t)\otimes\rho_{L}(t)\}$. The dissipator \mathcal{D}_{L}\big{(}\chi(t)\big{)} derived in Ref. Harbola et al. (2006) describes the action of the leads on the electronic system and has Lindblad form

[TABLE]

The Lindblad operators $s_{L}=|L\rangle\langle 0|$ and $s_{R}=|0\rangle\langle R|$ cause transitions between the states where an excess electron occupies either the left or right molecule of the dimer and the empty state of the system $|0\rangle$ where there are no excess charges. The Coulomb blockade is enforced by restricting the electronic subspace such that the identity is $I_{S}=|0\rangle\langle 0|+|L\rangle\langle L|+|R\rangle\langle R|$.

The hierarchical equations of motion can then be derived with respect to the system-bath coupling Tanimura (2006); Ishizaki and Fleming (2009); Kreisbeck et al. (2011); Tanimura (2012). A hierarchy of auxiliary density operators $\sigma^{(\mathbf{n})}(t)$ is formed which keep track of the system-bath correlations Zhu et al. (2012). The equations of motion are

[TABLE]

where the electronic system density matrix is given by $\rho_{S}(t)=\sigma^{(0)}(t)$. The auxiliary density matrices are labelled by the multi-index $\mathbf{n}=\{n_{L,0},...n_{L,K},n_{R,0},...n_{R,K}\}$ while coupling to higher and lower tiers is denoted by $\mathbf{n}\pm 1=\{n_{L,0},...n_{j,k}\pm 1,...,n_{R,K}\}$. The diagonal terms contain the Liouvillian \mathcal{L}\sigma^{(\mathbf{n})}(t)=-i[H_{S},\sigma^{(\mathbf{n})}(t)]+\mathcal{D}_{L}\big{(}\sigma^{(\mathbf{n})}(t)\big{)} which takes into account the unitary evolution of the auxiliary density matrices and the interaction with the leads. In Eq. (8), $j$ labels sites $|L\rangle$ and $|R\rangle$ coupled to the phonon bath while $k$ labels the Matsubara terms decaying with frequency $\nu_{k}$ in the expansion of the bath correlation function Shi et al. (2009). The coefficients of the correlation function $c_{j,k}$ depend on the form of the spectral density. In practice the hierarchy is truncated at a suitable level that ensures converged results while still enabling tractable numerical computation. Convergence can be improved through appropriate scaling of the auxiliary density matrices Shi et al. (2009) as well as through using a Markovian truncation term Ishizaki and Tanimura (2005) $\Xi=\sum_{j}\sum_{k=K}^{\infty}\frac{c_{jk}}{\nu_{k}}V^{\times}_{j}V^{\times}_{j}$ to reduce the number of Matsubara terms $K$ needed for convergence, thereby allowing investigation of quantum dynamics at low temperatures.

The hierarchical equations of motion can be formulated as a time-local master equation for the vector of auxiliary density matrices $\sigma(t)$

[TABLE]

with the hierarchy matrix $\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{J}$ acting as a time propagator. $\mathcal{H}_{0}$ describes the evolution of the system and auxiliary density matrices between electron jump events from the system to the drain lead while $\mathcal{H}_{J}$ describes these jump events and consists of the part of the dissipator \mathcal{D}_{L}\big{(}\rho_{S}(t)\big{)} containing $s_{R}$ in the top level of the hierarchy only.

We consider the current statistics between the system and the drain lead by augmenting $\mathcal{H}_{J}$ with a counting field $\chi$. The derivation of the counting statistics then follows in analogy with well known Markovian formalisms Bagrets and Nazarov (2003); Brandes (2005); Flindt et al. (2005); Marcos et al. (2010); Flindt et al. (2010). For ease of numerical computation we use the recursive scheme presented in Ref. Flindt et al. (2010). Introducing a perturbed time propagator $\mathcal{H}(\chi)=\mathcal{H}+\Delta\mathcal{H}(\chi)$ where $\Delta\mathcal{H}(\chi)=(e^{i\chi}-1)\mathcal{H}_{J}$, we obtain the $\chi$-dependent eigenvalue equation

[TABLE]

$|0(\chi)\rangle\rangle$ is the left steady state eigenvector of the hierarchy matrix and $\lambda_{0}(\chi)$ is the associated eigenvalue which tends to zero as $\chi\rightarrow 0$. The zero-frequency counting statistics are completely determined by this eigenvalue which can be expressed

[TABLE]

where $\langle\langle\tilde{0}|$ is the right steady state eigenvector of the hierarchy matrix.

Using the projectors $\mathcal{P}=|0\rangle\rangle\langle\langle\tilde{0}|$ and $\mathcal{Q}=I-\mathcal{P}$ the steady state can be written

[TABLE]

Additionally, using $\mathcal{H}=\mathcal{Q}\mathcal{H}\mathcal{Q}$ the eigenvalue equation becomes

[TABLE]

Introducing the pseudo-inverse $\mathcal{R}=\mathcal{Q}\mathcal{H}^{-1}\mathcal{Q}$ (more precisely this is the Drazin group generalized inverse Baiesi et al. (2009)) we obtain

[TABLE]

which on substituting into Eq. (12) gives us

[TABLE]

Expanding $\lambda_{0}(\chi)$, $|0(\chi)\rangle\rangle$ and $\Delta\mathcal{H}(\chi)$ about $\chi=0$ gives

[TABLE]

From this it is clear that the zero eigenvalue is equivalent to the cumulant generating function of the statistics. Substituting these expansions into Eqs. (11) and (15) allows calculation of the current cumulants to arbitrary order, though here we only make use of the first and second cumulants:

[TABLE]

[TABLE]

where the operators $\mathcal{H}^{(n)}=\mathcal{H}_{J}$.

To demonstrate the application of our non-perturbative framework we compare the electron counting statistics to a weak coupling approximation for the system-bath coupling Brandes (2005); Kießlich et al. (2007). The time generator for the perturbative approach is presented in Ref. Kießlich et al. (2007) and reads

[TABLE]

where

[TABLE]

are the general expressions for the dephasing rates with an arbitrary spectral density $J(\Delta)$ evaluated at the energy splitting between electronic eigenstates $\Delta=\sqrt{\epsilon^{2}+4T_{c}^{2}}$. Rates $\Gamma_{L}$ and $\Gamma_{R}$ represent the coupling to the left (source) and right (drain) leads while $\epsilon$ is the energy bias between left and right molecules and $T_{c}$ is the tunnelling coupling between the molecules. Initially we use a Drude-Lorentz spectral density $J(\omega)=\frac{2\lambda\Omega\omega}{\omega^{2}+\Omega^{2}}$ where $\lambda$ is the reorganisation energy of the bath quantifying the system-bath coupling strength, and $\Omega$ is the cutoff frequency quantifying the bath relaxation time. This form of the spectral density has been successfully used in combination with the hierarchical equations of motion to describe electron transfer in organic molecules Tanaka and Tanimura (2009); Wang et al. (2010). The central quantity we calculate for electron transport is the Fano factor $F=\frac{\langle\langle I^{2}\rangle\rangle}{\langle\langle I^{1}\rangle\rangle}$, which captures the deviation of the current fluctuations from Poissonian statistics.

The Fano factor as a function of dimer energy bias $\epsilon$ is shown in Fig. 2. For a double quantum dot system in the absence of a phonon bath the Fano factor is related to the steady state density matrix elements by the expression

[TABLE]

It is clear that the Fano factor is sensitive to both steady state populations and coherences with the contribution of the coherences dropping off as the energy bias becomes large. For small couplings to the drain lead ($\Gamma_{R}\ll 1$) this expression gives super-Poissonian current statistics ($F>1$) indicating charge transfer behaviour similar to the dynamical channel blockade Kießlich et al. (2007). In this regime delocalization of the charge across the dimer sets up two effective transport channels through either the symmetric or antisymmetric states of the dimer. A difference in the rate of transfer through each channel, in combination with the Coulomb blockade, leads to disordered transport and an increase in current fluctuations. This regime is sensitive to the effects of dephasing and so provides a good model in which to compare the effects of a phonon bath in both perturbative and non-perturbative theories. On introduction of the phonon bath, an asymmetry in the Fano factor with respect to positive or negative energy bias is obtained due to the increased rate of phonon emission relative to phonon absorption. However even for a small reorganisation energy ($\lambda\ll T_{c}$) there are clear differences between the predictions for the same set of parameters. In the weak system-bath coupling approximation, the interaction with the phonon bath induces dephasing processes which reduce the steady state coherence and consequently the Fano factor. Quantum dynamics in the non-perturbative regime is non-Markovian, where system-bath correlations help sustain electronic coherence for longer time scales. Consequently, a larger non-equilibrium steady state coherence is maintained, leading to larger Fano factors. The influence of this bath dynamics is most marked at higher temperatures ($\beta=0.1$), where a qualitative difference in the statistics for negative biases is observed. These results show that even in a regime where the weak-coupling approximation is expected to be valid to describe dynamics, the qualitative features of the zero-noise fluctuations may not be well-captured by these approximations.

The mean current and Fano factor as a function of reorganisation energy for low temperature ($\beta=0.4$) is shown in Fig. 3. For small reorganisations energies, where the bath has a weak influence on the electronic system, the weak coupling approximation and the non-perturbative theory give consistent results in both the first two current cumulants. As the reorganisation energy becomes the largest energy scale in the system, qualitative differences start to emerge in the predictions of both theories which signals the breakdown of the weak coupling approximation. This is due to the increased importance of system-bath correlations and the accompanying change in the nature of the electron-phonon relaxation process. The behaviour of both the mean and Fano factor is therefore more counter-intuitive, with the mean in particular exhibiting two maxima rather than the single maxima in the weak coupling approximation. This indicates the presence of two different optimization processes for large and small reorganisation energies with more moderate values leading to relaxation in a renormalized electronic basis such that population localization gives a reduced current output. The inset of Fig. 3 (b) demonstrates the convergence of the Fano factor with increasing Matsubara terms. As coupling to the bath increases more levels of the hierarchy and more Matsubara terms are required for convergence. At a certain reorganisation energy the computational effort will eventually become too great to achieve accurate results.

Coupling to discrete vibrational modes is also an important feature of electron transfer in molecular systems. Hence we investigate the effect of an underdamped Brownian oscillator spectral density on the current and Fano factor at high temperatures in Fig. 4. This spectral density has the form $J_{BO}(\omega)=\Theta(\omega)\frac{2\Lambda\gamma\omega_{0}^{2}\omega}{(\omega^{2}-\omega_{0}^{2})^{2}+\gamma^{2}\omega^{2}}$ where $\omega_{0}$ is the frequency of the mode, $\gamma$ is the damping constant and $\Lambda=S\omega_{0}$ is the reorganisation energy associated with the mode and $S$ is the Huang-Rhys factor. With this spectral density the hierarchical equations of motion take a slightly different form to that presented above Tanimura (2012). Rather than comparing the non-perturbative theory to the weak coupling approximation, a model in which a damped mode is coherently coupled to each site is used. In this coherent mode model the system Hamiltonian is

[TABLE]

where each mode is coupled to an individual site such that there is no direct interaction with the mode on the adjacent site. The coupling between the electronic and vibrational systems is given by $g_{j}=\sqrt{S_{j}}\omega_{j}$. Damping of the mode is provided using a Lindblad dissipator on the form

[TABLE]

where the quantities $\gamma_{j}$ are temperature independent damping rates and $n(\omega_{j})=(e^{\beta\omega_{j}}-1)^{-1}$ gives the occupation number of mode $j$ at inverse temperature $\beta$. Coupling with the leads is included in the same way in both models. The current and Fano factor are obtained using Eqs. 34 from Ref. Flindt et al. (2010) (the Markovian analogues of Eqs. (19) and (20) in this paper) with the same bath parameters used for both frameworks: $\omega_{0}=\omega_{L}=\omega_{R}$, $S=S_{L}=S_{R}$ and $\gamma=\gamma_{L}=\gamma_{R}$. Comparing the two models the current and Fano factor have the same qualitative features, with a peak at zero bias and enhancements in both quantities at resonance with the mode frequency. Asymmetry around zero bias is also present due to differences between phonon emission and absorption processes. The main qualitative difference however appears in the position of the peaks around the mode resonances. In the non-perturbative theory the peaks are shifted towards larger energies due to renormalisation of the electronic energy levels which is not fully captured in the coherent mode model.

A limitation to the non-perturbative and coherent mode models presented here lies in the validity of the description of the coupling between the electronic system and the leads. In the presence of strong electron-phonon interaction the leads may not interact with the electronic system independently of the phonon bath, rather with vibronic states. While this assumption has been used previously for counting statistics calculations in similar systems with arbitrary electron-phonon coupling Braggio et al. (2009b); Santamore et al. (2013), further developments to the non-perturbative theory will seek to investigate the extent of this issue.

In summary, we have have put forward a non-perturbative method to calculate electron counting statistics using the hierarchical equations of motion. Single molecule or single protein junctions in conjunction with single-electron counting devices Nishiguchi and Fujiwara (2009) will allow investigation of the statistics of the current passing through organic molecules thereby providing a powerful tool to complement transient spectroscopic studies. This in turn will lead to a deeper understanding of the microscopic mechanisms underlying electron transfer in a variety of organic systems while at the same time giving valuable insight for the design of reliable and better performing single-molecule devices.

Acknowledgements

The authors would like to thank Clive Emary for useful discussions. Financial support from the Engineering and Physical Sciences Research Council is gratefully acknowledged.