Generalized projected dynamics for non-system observables of non-equilibrium quantum impurity models

The study of open quantum impurity models, where the coupling of a small system to multiple baths drives it permanently away from the possibility of an equilibrium state, is an active and rapidly progressing field of research. It has recently become possible to make quantitative statements about experimentally measurable transport properties in certain cases [1–3]; however, in general many unresolved issues remain. For example, the nature of the charge and spin dynamics within the Kondo regime of a quantum dot driven out of equilibrium is currently under investigation [4], and basic questions regarding hysteresis and bi-stability in systems governed by strong electron–phonon couplings remain under debate [5–10]. Notably, many successful approaches to problems of this kind are based on master equation treatments and cumulant expansions [11–15], or on diagrammatic partial summations [16], all of which are approximate in general. A major theoretical challenge lies in the need to provide an accurate account of time propagation of open quantum systems, starting from some known initial state and proceeding all the way to an unknown steady-state.

Numerically exact methods play a particularly important role in the quest to obtain a reliable, unbiased description of non-equilibrium phenomena. Several different types of brute-force approaches developed in recent years have been applied to open non-equilibrium quantum systems. These include the time-dependent numerical renormalization group [17] and functional renormalization group [18–20], time-dependent density matrix renormalization group [21–24], iterative [25–28] and stochastic [29–33] diagrammatic methods and wavefunction based approaches [34, 35]. While the application of these approaches to the non-equilibrium Holstein, the Anderson impurity and the spin-fermion models has been very fruitful, they are still restricted to a relatively small range of parameters, typically characterized by a rapid decay to steady-state. Situations or observables exhibiting slow dynamics are inaccessible by these brute-force methods.

An alternative approach recently proposed by Cohen and Rabani [36] is based on a combination of a brute-force impurity solver (one of the above) with a generalized quantum master equation (GQME). The Nakajima–Zwanzig–Mori [37–39] formalism was used to derive an exact equation of motion for the reduced density matrix of the system, which includes a memory kernel giving rise to non-Markovian effects. This kernel, along with some information regarding the initial conditions, determines the dynamics of the system and contains all information about the time dependence of single-time system observables and their steady-state values. In many situations of interest, in particular when the bandwidth of the baths is large compared to other energy scales in the problem, the memory kernel is expected to decay rapidly to zero [4, 25, 36, 40]. Thus, one can safely truncate the memory kernel at a finite time, performing a 'cutoff approximation'. Brute-force impurity solvers limited to short times are well suited for the kernel's numerical evaluation up to the cutoff time, and once the memory kernel has been obtained, the GQME is exact and tractable at all times.

The GQME formalism has recently been combined with the Bold impurity solver [33, 41] to uncover the spin dynamics near the Kondo crossover temperature [4] and with the multilayer multi-configuration time-dependent Hartree method to reveal the nature of bi-stability in systems with electron–phonon couplings [40]. Despite the open nature of the systems studied in these works, transport properties were not addressed; this is due to one of the formalism's main limitations, in that observables outside the impurity part of the Hilbert space are not accessible, and only system observables such as the dot's population or magnetization can be calculated. On the other hand, perturbative expressions for transport properties in terms of vertex functions have been derived and evaluated before in approximate methodologies built on the GQME [42, 43], and it seems reasonable to expect that a general exact formulation in the spirit of [36] should exist.

In this paper we extend the GQME formalism (reviewed in section 2) to describe non-system observables. This allows for comparison of predictions made by the GQME with a much wider variety of experimental observables, of which an important example (worked out in detail here) is the current. Equally important, it facilitates access to the spectral functions by way of measuring the current to an infinitesimally coupled auxiliary bath [44–46]. The key idea discussed in section 3 is based on deriving a reduced equation of motion for the observable of interest, which can then be expressed in terms of the reduced density matrix. This leads to an introduction of an additional, observable-specific memory kernel with properties qualitatively similar to those of the memory kernel appearing in the standard GQME. Section 4 is devoted to expressing the steady-state properties in terms of the memory kernels alone, while in section 5 we show how the projected quantities appearing in the GQME can be translated into the language of ordinary observables expressed as second-quantized operators, using a non-interacting model as an illustrative example. In section 6 we present several test cases and examples for the non-interacting case, where the properties of memory kernels can be explored without the need for technically complicated numerical solvers. Finally, a summary is given in section 7.

We will begin by reviewing the derivation of exact projected (or reduced) equations of motion [47], and the process of going from projected to unprojected dynamics [48]. These details are provided here in a self-contained manner because they will be of particular importance later, when we discuss how the process can be generalized. Consider an operator Hilbert space $\mathcal {H}=\mathcal {S}\otimes \mathcal {B}$ composed of two subspaces $\mathcal {S}$ and $\mathcal {B}$, which we will call the system and bath subspaces. We are interested in a Hamiltonian of the form

$$\begin{equation} H=H_{\mathrm{S}}+H_{\mathrm{B}}+V, \end{equation} \tag{ 1 }$$

where $H_{\mathrm {S}}\in \mathcal {S}$ is the system or impurity Hamiltonian, $H_{\mathrm {B}}\in \mathcal {B}$ is the bath Hamiltonian and $V\in \mathcal {H},\, V\notin \mathcal {S},\,\mathcal {B},$ is the coupling Hamiltonian. Generally, the motivation for employing such a description is to describe a small, strongly interacting impurity coupled to large non-interacting baths, but we need make no further assumptions at this stage. We can now define a projection operator P onto the $\mathcal {S}$ subspace by tracing out the bath degrees of freedom, in the process also defining its complementary operator Q:

$$\begin{eqnarray} P & = & \rho_{\mathrm{B}}\mathrm{Tr}_{\mathrm{B}}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} Q & = & 1-P. \end{eqnarray} \tag{ 3 }$$

Here ρ_B = e^−βH_B/Tr_B e^−βH_B. We also define $\rho _{\mathrm {S}}\in \mathcal {S}$ to be the initial impurity density matrix, and ρ₀ = ρ_B⊗ρ_S the initial full density matrix. The expectation value of a system operator $A\in \mathcal {S}$ is given by

$$\begin{eqnarray} \left\langle A(t)\right\rangle & = & \mathrm{Tr}\rho(t)A \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\hphantom{\left\langle A(t)\right\rangle} = \mathrm{Tr}_{\mathrm{S}}\left[\left(\mathrm{Tr}_{\mathrm{B}}\rho(t)\right)A\right] \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\hphantom{\left\langle A(t)\right\rangle}\equiv \mathrm{Tr}_{\mathrm{S}}\left\{ \sigma(t)A\right\} \end{eqnarray} \tag{ 6 }$$

and the reduced density matrix σ(t) = Tr_Bρ(t) contains information about all single-time properties of system observables. This object has a lower dimensionality than that of ρ, and it would thus be economical to describe its equations of motion without referring to the system as a whole. This is the basic idea behind reduced quantum dynamics.

To proceed, one considers the Liouville–von Neumann equation, which governs the dynamics of the full density matrix:

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}\rho & =\left[H,\rho\right]\equiv & \mathcal{L}\rho. \end{eqnarray} \tag{ 7 }$$

The Liouvillian superoperator $\mathcal {L}$ denotes performing a commutation with the Hamiltonian, such that $\mathcal {L}A\equiv \left [H,A\right ]$. We also define $\mathcal {L}_{\mathrm {S}}A\equiv \left [H_{\mathrm {S}},A\right ]$, $\mathcal {L}_{V}A\equiv \left [V,A\right ]$ and $\mathcal {L}_{\mathrm {B}}A\equiv \left [H_{\mathrm {B}},A\right ]$. Applying each of the projection operators from the left and using 1 = P + Q within the commutator gives

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}P\rho & = & P\left[H,\left(P+Q\right)\rho\right], \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}Q\rho & = & Q\left[H,\left(P+Q\right)\rho\right]. \end{eqnarray} \tag{ 9 }$$

Equation (9) has the formal solution

$$\begin{eqnarray} &&Q\rho = \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\rho_{0} -\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\, \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}\tau}Q\mathcal{L}\rho_{\mathrm{B}}\sigma\left(t-\tau\right), \end{eqnarray} \tag{ 10 }$$

which can be inserted into equation (8) to obtain the Nakajima–Zwanzig–Mori equation (NZME) [37–39]:

$$\begin{eqnarray} \mathrm{i}\hbar\dot{\sigma}(t) & = & \mathcal{L}_{\mathrm{S}}\sigma(t)+\vartheta(t)-\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\kappa(\tau)\sigma\left(t-\tau\right), \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \kappa(t) & \equiv & \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\mathcal{L}\rho_{\mathrm{B}}\right\} , \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \vartheta(t) & \equiv & \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\rho_{0}\right\} . \end{eqnarray} \tag{ 13 }$$

Let us take a moment to examine the important relation equation (11). It has the form of an operator linear Volterra integro-differential equation of the second kind. As we have made no approximations, it is exact; yet it contains only operators and superoperators within the low-dimensional system space. The time derivative of the reduced density matrix σ is given by the sum of three contributions: the first term ($\mathcal {L}_{\mathrm {S}}\sigma (t)$) describes the exact evolution of the system if the coupling to the bath were set to zero. The second term (ϑ(t)) expresses initial correlations between the system and bath, and it is easy to verify from its definition in equation (13) that it equals zero for the factorized initial conditions ρ₀ = ρ_B⊗ρ_S (we will assume this later, but keep this term for generality, as access to general initial conditions is of some interest when considering, for instance, quenching). The last term includes the memory kernel ($\kappa (\tau )$), and depends on the complete history of σ(t) at earlier times. The appearance of this non-Markovian term is the price of going to reduced dynamics, and to make headway with the NZME one must begin by evaluating κ(t).

The definition of κ(t) in equation (12) includes the troublesome component $\mathrm {e}^{-\frac {\mathrm {i}}{\hbar }Q\mathcal {L}\tau }$. To understand why it is troubling, consider the following: it is easy to show that the superoperator $\mathrm {e}^{-\frac {\mathrm {i}}{\hbar }\mathcal {L}\tau }$ evolves the density matrix with respect to the Hamiltonian, thus simply expressing our familiar notion of dynamics:

$$\begin{equation} \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}i\mathcal{L}\tau}\rho=\mathrm{e}^{\frac{\mathrm{i}}{\hbar}H\tau}\rho \,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}H\tau}. \end{equation} \tag{ 14 }$$

The modified operator $\mathrm {e}^{-\frac {\mathrm {i}}{\hbar }Q\mathcal {L}\tau }$, however, contains a projection operator within the exponent, and so does something else entirely—something which turns out to be substantially harder to understand or calculate. Our next step is therefore to get rid of these inconvenient projected dynamics. While several ways to go about this task exist, we will limit the discussion to a particular method suggested by Zhang et al [48].

Consider the function ϑ(t) of equation (13). By applying the identity

$$\begin{equation} \mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}=\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}+\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau \,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}\left(t-\tau\right)}P\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}\tau} \end{equation} \tag{ 15 }$$

to its definition, we can obtain

$$\begin{eqnarray} \vartheta(t) & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\rho_{0}\right\} \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\hphantom{\vartheta(t)} = \mathrm{Tr}_{\mathrm{B}}\Bigg\{\mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}Q\rho_{0} +\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}\left(t-\tau\right)}P\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}\tau}Q\rho_{0}\Bigg\} \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&\hphantom{\vartheta(t)}= \Xi(t)-\Phi(t)\sigma(0) +\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\Phi\left(t-\tau\right)\vartheta(\tau), \end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray} \Xi(t) & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}\rho_{0}\right\} , \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \Phi(t) & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}\rho_{\mathrm{B}}\right\} . \end{eqnarray} \tag{ 20 }$$

Applying the same identity equation (15) to equation (12) yields

$$\begin{equation} \kappa(t)=\mathrm{i}\hbar\dot{\Phi}(t)-\Phi(\tau)\mathcal{L}_{\mathrm{S}}+\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\Phi\left(t-\tau\right)\kappa(\tau). \end{equation} \tag{ 21 }$$

Equations (18) and (21) are superoperator linear Volterra integral equations of the second kind, with both the inhomogeneous contributions and the kernels determined by the combination of equations (19) and (20) and the form of the system Liouvillian operator. Like the NZME, equation (11), they consist of objects which inhabit the low-dimensional impurity subspace—however, they have a higher dimensionality due to their superoperator nature (if σ can be represented by an N × N matrix, then ϑ and κ are N² × N²). Their importance lies in the fact that Φ and Ξ, which are propagated by the full Hamiltonian with normal dynamics, can be written in terms of physical observables; this means they can be evaluated with a variety of computational methods, and then used to solve equations (18) and (21) numerically.

We now have the necessary machinery at hand to introduce the cutoff approximation: if we have some way of evaluating κ(t) up to some finite time, it is sometimes possible to make an ansatz about later times. Importantly, if the memory has decayed to zero to within a numerical accuracy over a finite time, one can assume that it will remain zero at all later times. One then solves the NZME with this cutoff memory kernel to obtain an approximate value for σ(t); however, if σ(t) can be converged in the cutoff time to within the desired accuracy, the entire procedure is numerically exact. Note that while in principle this procedure can be performed for any Hamiltonian (regardless of the form of the interactions), for it to be beneficial in practice the system in question should exhibit dynamical timescales substantially longer than those of the memory decay time.

The time-dependent electronic current flowing through an impurity does not have a single definition, as it depends in general on the topology of the surface through which electronic flow is measured. In steady state populations must be constant, and the current must therefore become independent of this definition (if it is unique); however, the definition itself remains arbitrary. This is well known and usually does not warrant much discussion, yet in the context of reduced dynamics a subtle point occurs: if the impurity model in question is, for instance, given by a chain Hamiltonian, current may be measured at any point along the chain and may be obtained from knowledge of σ(t). However, in models where not all current must flow between impurity sites, it is necessary to measure currents at the junction between the impurity and one of the leads (baths). In a setup involving two Fermionic baths held at different chemical potentials, often referred to as the left (L) and right (R) leads, one is therefore interested in quantities such as the so-called 'left current':

$$\begin{equation} I_{L}\equiv\frac{\mathrm{d}}{\mathrm{d}t}eN_{L}=\sum_{q\in L}\left\langle \frac{\mathrm{i}e}{\hbar}\left[H,a_{q}^{\dagger}a_{q}\right]\right\rangle . \end{equation} \tag{ 22 }$$

Here the a^†_k and a_k are creation and destruction operators in the left lead subspace $L\subset \mathcal {B}$. The current operator will therefore in general not be a member of the impurity subspace $\mathcal {S}$, and cannot be obtained from equation (6) along with knowledge of the reduced density matrix σ(t). It should be mentioned briefly that one simple way of dealing with this issue is to define an effective Hamiltonian and a repartitioning of $\mathcal {H}$ in such a way that the current can be measured within the system; however, this will invariably raise the dimensionality of $\mathcal {S}$, which may complicate the problem beyond the applicability of many numerical methods.

In order to allow for the calculation of non-system observables, we proceed by deriving a Nakajima–Zwanzig–Mori-like equation for the expectation value of a general operator I. While I implies that we are interested in the current, nothing in this section is limited to that specific case. The ideas to follow could work equally well for any operator, but for the current one might expect them to converge quickly with a cutoff time, since it is expected to be determined largely by quantities local to the dot.

We start from the equation of motion for I in the Schrödinger picture, where ρ has a time dependence but I does not:

$$\begin{equation} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}I\rho=I\left[H,\rho\right]. \end{equation} \tag{ 23 }$$

Applying the projection operators using the definitions and procedure of the previous section yields

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}PI\rho & = & PI\mathcal{L}P\rho+PI\mathcal{L}Q\rho, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}QI\rho & = & QI\mathcal{L}P\rho+QI\mathcal{L}Q\rho. \end{eqnarray} \tag{ 25 }$$

In addition to these, we still have equations (8) and (9) for the density matrix, which may be written in the form

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}P\rho & = & P\mathcal{L}P\rho+P\mathcal{L}Q\rho, \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}Q\rho & = & Q\mathcal{L}P\rho+Q\mathcal{L}Q\rho. \end{eqnarray} \tag{ 27 }$$

The latter equation, as before, has the formal solution equation (10). Putting this expression together with equation (24) and defining ι(t) = Tr_BIρ, we obtain

$$\begin{eqnarray} &&\fl \mathrm{i}\hbar \frac{{\rm d}\iota(t)}{{\rm d}t} = \mathrm{Tr}_{\mathrm{B}}I\mathcal{L}\rho_{\mathrm{B}}\sigma(t)+\mathrm{Tr}_{\mathrm{B}}I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\rho_{0} -\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\,\mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}\left(t-\tau\right)}Q\mathcal{L}\rho_{\mathrm{B}}\right\} \sigma(\tau) \end{eqnarray} \tag{ 28 }$$

or

$$\begin{equation} \mathrm{i}\hbar\frac{{\rm d}\iota(t)}{{\rm d}t}=\mathcal{L}_{\iota}\sigma(t)+\vartheta_{\iota}(t)-\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\,\kappa_{\iota}\left(t-\tau\right)\sigma(\tau). \end{equation} \tag{ 29 }$$

The terms of this equation appear similar to those of the NZME in equation (11), though it is a solution in closed form rather than an integro-differential equation, since ι(t) appears only on the left-hand side. In writing it we have defined

$$\begin{eqnarray} \mathcal{L}_{\iota} & \equiv & \mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\rho_{\mathrm{B}}\right\} , \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \vartheta_{\iota}(t) & \equiv & \mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\rho_{0}\right\} , \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} \kappa_{\iota}(t) & \equiv & \mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}t}Q\mathcal{L}\rho_{\mathrm{B}}\right\} . \end{eqnarray} \tag{ 32 }$$

The initial correlation term ϑ_ι is once again zero for uncorrelated initial conditions and this time we will remove it for the sake of brevity. As in the formalism for σ, in order to phrase everything in terms of quantities with unprojected dynamics we now once again perform the Zhang–Ka–Geva transformation [48] on the current memory kernel κ_ι. Applying the identity (15) allows us to write

$$\begin{eqnarray} &&\fl \kappa_{\iota}(t) = \mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}Q\mathcal{L}\rho_{\mathrm{B}}\right\} +\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}\left(t-\tau\right)}P\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Q\mathcal{L}\tau}Q\mathcal{L}\rho_{\mathrm{B}}\right\} \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} &&\fl \hphantom{\kappa_{\iota}(t) } = \mathrm{i}\hbar\dot{\Phi}_{\iota}(t)-\Phi_{\iota}(t)\mathcal{L}_{\mathrm{S}} +\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\Phi_{\iota}\left(t-\tau\right)\kappa(\tau) \end{eqnarray} \tag{ 34 }$$

or

$$\begin{eqnarray} \kappa_{\iota}(t) & = & \mathrm{i}\hbar\dot{\Phi}_{\iota}(t)-\Phi_{\iota}(t)\mathcal{L}_{\mathrm{S}} +\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\Phi_{\iota}\left(t-\tau\right)\kappa(\tau). \end{eqnarray} \tag{ 35 }$$

Once again, this is a closed form solution rather than an integral equation. Its inputs are the same $\kappa (\tau )$ defined in equation (12), as well as the new quantity

$$\begin{eqnarray} \Phi_{\iota}(t) & = & \mathrm{Tr}_{\mathrm{B}}\left\{ I\mathcal{L}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}\rho_{\mathrm{B}}\right\} . \end{eqnarray} \tag{ 36 }$$

Equations (35) and (29) amount to a generalization of the Nakajima–Zwanzig–Mori formalism to non-system operators. The structure of these equations is reminiscent of the structure of the corresponding equations in the original theory, on which the extension relies—and yet they are simpler in a certain sense, as they are closed form solutions up to quadrature rather than integro-differential or integral equations. In addition to σ(t) and κ(t), which can be obtained from the original theory, the extended formalism relies on a new input, Φ_ι(t), which is defined in terms of regular (rather than projected) time propagation and must be calculated explicitly. Once Φ_ι(t) is available one can solve equation (35) to obtain κ_ι(t), and then solve equation (29) to obtain ι(t), a system-space operator which can be traced over to obtain the expectation value of the operator I.

If we wish to examine the t → ∞ limit of σ(t), it is more convenient to define the Laplace transform

$$\begin{equation} \hat{\sigma}\left(z\right)=\int_{0}^{\infty}\mathrm{e}^{-zt}\sigma(t)\,\mathrm{d}t. \end{equation} \tag{ 37 }$$

When applied to equation (11), this yields

$$\begin{eqnarray} \mathrm{i}\hbar\left[z\hat{\sigma}\left(z\right)-\sigma(0)\right] & = & \mathcal{L}_{\mathrm{S}}\hat{\sigma}\left(z\right)+\hat{\vartheta}\left(z\right)-\frac{\mathrm{i}}{\hbar}\hat{\kappa}\left(z\right)\hat{\sigma}\left(z\right), \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&\Downarrow\nonumber \\ &&\hat{\sigma}\left(z\right) = \frac{\sigma(0)+\frac{1}{\mathrm{i}\hbar}\hat{\vartheta}\left(z\right)}{z-\frac{1}{\mathrm{i}\hbar}\mathcal{L}_{\mathrm{S}}+\frac{1}{\hbar^{2}}\hat{\kappa}\left(z\right)}. \end{eqnarray} \tag{ 39 }$$

Using the final value theorem $\sigma \left (\infty \right )=\lim _{z\rightarrow 0}z\hat {\sigma }\left (z\right )$, we can obtain an expression for $\hat {\sigma }$ at long times:

$$\begin{equation} \sigma\left(t\rightarrow\infty\right)=\lim_{z\rightarrow0}\frac{\mathrm{i}\hbar\sigma(0)+\hat{\vartheta}\left(z\right)}{\left(\mathrm{i}\hbar+\frac{\mathrm{i}}{\hbar}\frac{1}{z}\hat{\kappa}\left(z\right)-\frac{1}{z}\mathcal{L}_{\mathrm{S}}\right)}. \end{equation} \tag{ 40 }$$

We can also obtain a stationary-state equation by considering a time independent solution $\sigma \left (t\rightarrow \infty \right )$ to equation (11), such that we can set the time derivative to zero and take σ outside the integral before taking the Laplace transform. If we also assume that the initial correlations are either zero to begin with or die out at infinite time, this gives

$$\begin{equation} \left(\mathcal{L}_{\mathrm{S}}-\frac{\mathrm{i}}{\hbar}\hat{\kappa}\left(z\rightarrow \mathrm{i}0\right)\right)\sigma\left(t\rightarrow\infty\right)=0. \end{equation} \tag{ 41 }$$

This last equation is of particular interest, because it allows us to go from the memory kernel and system Liouvillian directly to the steady state properties of the reduced density matrix, without passing through the dynamics and without any reference to the initial state or correlations of the system. This is very useful when we are interested in general questions regarding the steady state, such as that of its existence or uniqueness. When applying the cutoff approximation, $\hat {\kappa }\left (z\rightarrow \mathrm {i}0\right )$ must be calculated to sufficient accuracy that the steady state density matrix converges.

It is natural to attempt deriving a similar expression for the current directly at steady state. One way of going about this task is to begin with equation (29) and take the Laplace transform:

$$\begin{eqnarray} \mathrm{i}\hbar\left[z\hat{\iota}\left(z\right)-\overset{{\scriptstyle =0}}{\overbrace{\iota(0)}}\right] & = & \mathcal{L}_{\iota}\hat{\sigma}\left(z\right)+\hat{\vartheta}_{\iota}\left(z\right) -\frac{\mathrm{i}}{\hbar}\hat{\kappa}_{\iota}\left(z\right)\hat{\sigma}\left(z\right). \end{eqnarray} \tag{ 42 }$$

Extracting $z\hat {\iota }\left (z\right )$ and using the final value theorem then gives

$$\begin{eqnarray} &&\iota\left(t\rightarrow\infty\right) = \lim_{z\rightarrow0}\frac{1}{\mathrm{i}\hbar}\left(\mathcal{L}_{\iota}-\frac{\mathrm{i}}{\hbar}\hat{\kappa}_{\iota}\left(z\right)\right)\hat{\sigma}\left(z\right) \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} &&\hphantom{\iota\left(t\rightarrow\infty\right)} = \lim_{z\rightarrow0}\frac{1}{\mathrm{i}\hbar z}\left(\mathcal{L}_{\iota}-\frac{\mathrm{i}}{\hbar}\hat{\kappa}_{\iota}\left(z\right)\right)\sigma\left(t\rightarrow\infty\right). \end{eqnarray} \tag{ 44 }$$

This suggests that in order for a steady state to exist, we must have

$$\begin{equation} \lim_{z\rightarrow0}\mathcal{L}_{\iota}-\frac{\mathrm{i}}{\hbar}\hat{\kappa}_{\iota}\left(z\right)\sim z. \end{equation} \tag{ 45 }$$

The constant of proportionality (itself a superoperator) determines the value of the current at steady state. In the case of impurity observables, it is sufficient to know the zero frequency component of the memory kernel in order to obtain the steady-state value. Here, however, one must also know something about the low-frequency properties (or linear frequency response) of the current memory kernel, $\hat {\kappa }_{\iota }\left (z\right )$.

Everything up to this point has been independent of the details of any particular model, requiring only that a partitioning between the impurity and bath part be made. In order to illustrate the process of using the formalism presented above in a particular model, we will continue by way of the simplest possible example: that of a non-interacting junction (the formalism is not limited to this case [4, 36, 40]). This model, often called the resonant level model, is defined by the Hamiltonian

$$\begin{eqnarray} H & = & H_{\mathrm{S}}+H_{L}+V, \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} H_{\mathrm{S}} & = & \varepsilon d^{\dagger}d, \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} H_{L} & = & \sum_{q}\varepsilon_{q}a_{q}^{\dagger}a_{q}, \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} V & = & \sum_{q}t_{q}da_{q}^{\dagger}+t_{q}^{*}a_{q}d^{\dagger}. \end{eqnarray} \tag{ 49 }$$

A complete definition must include the ε_q and t_q, and in this case all necessary information is contained in the lead coupling function

$$\begin{equation} \Gamma\left(\omega\right)=2\pi\sum_{q}\left|t_{q}\right|^{2}\delta\left(\omega-\omega_{q}\right). \end{equation} \tag{ 50 }$$

The first step in the calculation is the evaluation of the system Liovillian. It is convenient to work in the Hubbard representation for operators in the impurity subspace: an operator $\skew3\hat {A}\in \mathcal {S}$ can be written as $\hat {A}=\sum _{ij}a_{ij}\left |i\right \rangle \left \langle j\right |$, where the indices i and j can take on the values of states in the impurity subspace—in this case 0 and 1 for unoccupied and occupied, respectively. This superoperator simply performs a commutation with the system Hamiltonian, and using equation (47) to insert the explicit form of H_S yields an expressions for $\mathcal {L}_{\mathrm {S}}$ in matrix (or tetradic) form

$$\begin{eqnarray} \left[\mathcal{L}_{\mathrm{s}}\right]_{ij,kl} & = & \mathrm{Tr}_{\mathrm{S}}\left\{ \left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}\mathcal{L}_{\mathrm{S}}\left|k\right\rangle \left\langle l\right|\right\} \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} &&\hphantom{\left[\mathcal{L}_{\mathrm{s}}\right]_{ij,kl} } = \sum_{m=0}^{1}\left\langle m\right|\left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}\left[\varepsilon\left|1\right\rangle \left\langle 1\right|,\left|k\right\rangle \left\langle l\right|\right]\left|m\right\rangle \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} &&\hphantom{\left[\mathcal{L}_{\mathrm{s}}\right]_{ij,kl} } = \varepsilon\left[\delta_{jl0}\delta_{ik1}+\delta_{ijkl1}-\delta_{jl1}\delta_{ik}\right]. \end{eqnarray} \tag{ 53 }$$

Here δ_a1a2...aN is one if all indices take the same value and zero otherwise. When no bath is present the model is reduced to a two-level system, and equation (11) gives the expected result

$$\begin{eqnarray} \mathrm{i}\hbar\frac{\mathrm{d}\sigma_{ij}}{\mathrm{d}t} & = & \sum_{kl}\left[\mathcal{L}_{\mathrm{s}}\right]_{ij,kl}\sigma_{kl} \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} &&\hphantom{\mathrm{i}\hbar\frac{\mathrm{d}\sigma_{ij}}{\mathrm{d}t} } = \varepsilon\left(\delta_{i1}\delta_{j0}\sigma_{10}-\delta_{i0}\delta_{j1}\sigma_{01}\right). \end{eqnarray} \tag{ 55 }$$

That is, off-diagonal density matrix elements oscillate with a frequency $\frac {\varepsilon }{\hbar }$ while diagonal elements remain stationary.

Next, we need to evaluate the memory kernel. This requires the evaluation of the superoperator

$$\begin{eqnarray} &&\Phi(t)A = \mathrm{Tr}_{\mathrm{B}}\left\{ \mathcal{L}_{V}\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}\mathcal{L}t}\rho_{\mathrm{B}}A\right\} \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} \hphantom{\Phi(t)A}& = & \mathrm{Tr}_{\mathrm{B}}\left\{ V\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Ht}\rho_{\mathrm{B}}A\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar}Ht} -\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Ht}\rho_{\mathrm{B}}A\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar}Ht}V\right\} \end{eqnarray} \tag{ 57 }$$

which can also be represented in matrix form

$$\begin{eqnarray} \Phi_{ij,kl}(t) & = & \mathrm{Tr}_{\mathrm{S}}\left\{ \left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}\phi(t)\left|k\right\rangle \left\langle l\right|\right\} \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} &&\hphantom{\Phi_{ij,kl}(t)}= A-A^{\prime} \end{eqnarray} \tag{ 59 }$$

with

$$\begin{eqnarray} A & \equiv & \mathrm{Tr}_{\mathrm{S}}\left\{ \left|j\right\rangle \left\langle i\right|\mathrm{Tr}_{\mathrm{B}}\left\{ V\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Ht}\rho_{\mathrm{B}}\left|k\right\rangle \left\langle l\right|\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar}Ht}\right\} \right\} , \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} A^{\prime} & \equiv & \mathrm{Tr}_{\mathrm{S}}\left\{ \left|j\right\rangle \left\langle i\right|\mathrm{Tr}_{\mathrm{B}}\left\{\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Ht}\rho_{\mathrm{B}}\left|k\right\rangle \left\langle l\right|\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar}Ht}V\right\} \right\} . \end{eqnarray} \tag{ 61 }$$

Consider the term A. Let us perform the trace over the impurity space and take the Hubbard operators to second quantized form

$$\begin{eqnarray} &&\fl A = \mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\mathrm{e}^{\frac{\mathrm{i}}{\hbar}Ht}\right. \,\left(\delta_{i1}\delta_{j1}d^{\dagger}d+\delta_{i0}\delta_{j0}dd^{\dagger}+\delta_{i1}\delta_{j0}d+\delta_{i0}\delta_{j1}d^{\dagger}\right)\left.V\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar}Ht}\left|k\right\rangle \right\} \end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray} &&\fl \hphantom{A} = \mathrm{Tr}_{\mathrm{B}}\Big\{\rho_{\mathrm{B}}\left\langle l\right| \,\left(\delta_{i1}\delta_{j1}d^{\dagger}(t)d(t)+\delta_{i0}\delta_{j0}d(t)d^{\dagger}(t)\right. \,\left.+\delta_{i1}\delta_{j0}d(t)+\delta_{i0}\delta_{j1}d^{\dagger}(t)\right) \, V(t)\left|k\right\rangle \Big\}. \end{eqnarray} \tag{ 63 }$$

In the final step, the operators were given their full time dependence in the Heisenberg picture. Using $V(t)=\sum _{q}t_{q}d(t)a_{q}^{\dagger }(t)+t_{q}^{*}a_{q}(t)d^{\dagger }(t)$ and the fact that all pairs of dot and lead operators maintain normal commutation relations when taken at the same times, one can now show that

$$\begin{eqnarray} &&\fl A = \sum_{q}t_{q}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i0}\delta_{j0}d(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} +\sum_{q}t_{q}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i0}\delta_{j1}d^{\dagger}(t)d(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} \nonumber \\ &&\fl\hphantom{A =} -\sum_{q}t_{q}^{*}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i1}\delta_{j1}d^{\dagger}(t)a_{q}(t)\left|k\right\rangle \right\}-\sum_{q}t_{q}^{*}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i1}\delta_{j0}d(t)d^{\dagger}(t)a_{q}(t)\left|k\right\rangle \right\} .\nonumber\\ \end{eqnarray} \tag{ 64 }$$

Similarly,

$$\begin{eqnarray} &&\fl A^{\prime} = \sum_{q}t_{q}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i1}\delta_{j1}d(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} -\sum_{q}t_{q}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i0}\delta_{j1}d(t)d^{\dagger}(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} \nonumber \\ &&\fl\hphantom{A' =}-\sum_{q}t_{q}^{*}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i0}\delta_{j0}d^{\dagger}(t)a_{q}(t)\left|k\right\rangle \right\} +\sum_{q}t_{q}^{*}\mathrm{Tr}_{\mathrm{B}}\left\{ \rho_{\mathrm{B}}\left\langle l\right|\delta_{i1}\delta_{j0}d^{\dagger}(t)d(t)a_{q}(t)\left|k\right\rangle \right\}. \nonumber \\ \end{eqnarray} \tag{ 65 }$$

Putting the expressions for A and A' into their defining equation then yields

$$\begin{eqnarray} \Phi_{ij,kl}(t) & = & -2\mathrm{i}\left(\delta_{i1}\delta_{j1}-\delta_{i0}\delta_{j0}\right)\Im\left\{ \varphi_{kl}\right\} +\delta_{i0}\delta_{j1}\psi_{kl}-\delta_{i1}\delta_{j0}\psi_{lk}^{*} \end{eqnarray} \tag{ 66 }$$

with

$$\begin{eqnarray} \varphi_{kl} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q}t_{q}\rho_{\mathrm{B}}\left\langle l\right|d(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} , \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} \psi_{kl} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q}t_{q}\rho_{\mathrm{B}}\left\langle l\right|a_{q}^{\dagger}(t)\left|k\right\rangle \right\} . \end{eqnarray} \tag{ 68 }$$

The φ elements have a rather simple physical interpretation: they are directly proportional to the time derivative of the total population on the dot. The ψ elements are harder to interpret in such a manner.

The equations therefore collapse to a simple form, phrased in terms of the system-space matrix elements φ_kl and ψ_kl of normal second quantization operators propagated under the influence of the full Hamiltonian. Some further simplification can be made by considering equation (67) as a function of time when we go to the interaction picture under H₀ = H_S + H_B:

$$\begin{eqnarray} \varphi_{kl} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q}t_{q}\rho_{\mathrm{B}}\left\langle l\right|\mathrm{e}^{\mathrm{i}Ht}da_{q}^{\dagger}\,\mathrm{e}^{-\mathrm{i}Ht}\left|k\right\rangle \right\} \end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray} &&\hphantom{\varphi_{kl} }= \sum_{q}t_{q}\mathrm{Tr}_{\mathrm{B}} \,\left\{ \rho_{\mathrm{B}}\left\langle l\right|U^{\dagger}(t)d_{H_{0}}(t)a_{H_{0},q}^{\dagger}(t)U(t)\left|k\right\rangle \right\} . \end{eqnarray} \tag{ 70 }$$

It is easy to verify that the time dependence of d and a^†_q in the interaction-picture is described by a simple oscillation, and that the U and U^† are sums over products of terms containing either a^†_qd or d^†a_q in the interaction picture. The trace over the bath may then be performed at time zero, throwing out all terms which do not have the same number of a_q and a^†_q operators. Yet, from the argument we have just stated, such terms will also have the same number of d and d^† operators, and φ_kl must be zero unless k = l. For similar considerations, ψ_kl is zero unless k ≠ l and we have

$$\begin{eqnarray} \varphi_{kl} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q}t_{p}\rho_{\mathrm{B}}\left\langle l\right|d(t)a_{q}^{\dagger}(t)\left|k\right\rangle \right\} \delta_{kl}, \end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray} \psi_{kl} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q}t_{q}\rho_{\mathrm{B}}\left\langle l\right|a_{q}^{\dagger}(t)\left|k\right\rangle \right\} \left(1-\delta_{kl}\right). \end{eqnarray} \tag{ 72 }$$

Considering the role of φ and ψ in equation (66), one can see that Φ contains terms which couple the populations and terms which couple the coherences; it contains no terms which couple the populations to the coherences. Examining equation (21), one realizes that terms with no inhomogeneous contribution must identically vanish, such that

$$\begin{equation} \kappa_{ij,kl}=\left\{\begin{array}{@{}ll@{}} \kappa_{ii,kk} & i=j,\, k=l,\\ \kappa_{ij,kl} & i\neq j,\, k\neq l,\\ 0 & \mathrm{otherwise}. \end{array}\right. \end{equation} \tag{ 73 }$$

Similarly, equation (11), along with the Liouvillian (53), immediately leads us to the conclusion that the diagonal elements of σ form one coupled block within the formalism, while the off-diagonal elements of σ form a second block: in other words, within the resonant level model, the reduced dynamics of the diagonal elements (the populations) are decoupled from those of the off-diagonal elements (the coherences).

Among other things, this implies that if we are interested only in the populations we do not need to calculate the ψ_kl, and vice-versa for the coherences. Since the populations are also unaffected by the Liouvillian, and assuming factorized initial conditions, the equation of motion turns from a superoperator equation into a matrix equation for the population vector σ_ii:

$$\begin{equation} \mathrm{i}\hbar\frac{\mathrm{d}}{\mathrm{d}t}\sigma_{ii}(t)=-\frac{\mathrm{i}}{\hbar}\int_{0}^{t}\mathrm{d}\tau\,\sum_{j}\kappa_{ii,jj}(\tau)\sigma_{jj}\left(t-\tau\right). \end{equation} \tag{ 74 }$$

Interestingly, this analytic block structure conclusion holds in the generalized Holstein model [40] as well (though this will not be shown here), and continues to hold for the Anderson model [36] even in the presence of a magnetic field [4].

We continue to examine the memory formalism for the non-system current operator in the non-interacting case. We will define a simplified left current operator as

$$\begin{equation} \tilde{I}\equiv\sum_{q\in L}t_{q}da_{q}^{\dagger} \end{equation} \tag{ 75 }$$

such that $\tilde {\iota }(t)=\mathrm {Tr}_{\mathrm {B}}\left \{ \tilde {I}\rho (t)\right \}$ and the physical current is

$$\begin{equation} \left\langle I(t)\right\rangle =-2\frac{e}{\hbar}\Im\left\langle \tilde{I}(t)\right\rangle =-2\frac{e}{\hbar}\Im\mathrm{Tr}_{\mathrm{S}}\left\{ \tilde{\iota}(t)\right\} . \end{equation} \tag{ 76 }$$

First, consider the current Liouvillian operator:

$$\begin{eqnarray} \mathrm{Tr}_{\mathrm{B}}\left\{ \tilde{I}\mathcal{L}\rho_{\mathrm{B}}\right\} & = & \mathrm{Tr}_{\mathrm{B}}\left\{ \tilde{I}\mathcal{L}\rho_{\mathrm{B}}\right\} \end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray} &&\hphantom{\mathrm{Tr}_{\mathrm{B}}\left\{ \tilde{I}\mathcal{L}\rho_{\mathrm{B}}\right\} }= \mathrm{Tr}_{\mathrm{B}}\left\{ \sum_{q\in L}t_{q}\,\mathrm{d}a_{q}^{\dagger}\mathcal{L}_{V}\rho_{\mathrm{B}}\right\} . \end{eqnarray} \tag{ 78 }$$

The term containing $\mathcal {L}_{\mathrm {B}}$ can be dropped because $\mathcal {L}_{\mathrm {B}}\rho _{\mathrm {B}}A=\left (\mathcal {L}_{\mathrm {B}}\rho _{\mathrm {B}}\right )A=0$ for any system variable A, while the term containing $\mathcal {L}_{\mathrm {S}}$ is zero because it must contain a trace over an odd number of lead creation or destruction operators. It is then simple to show by the same procedure from which the system Liouvillian was derived that

$$\begin{eqnarray} &&\fl \left[\mathcal{L}_{\iota}\right]_{ij,ml} = \mathrm{Tr}_{\mathrm{S}}\left\{ \left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}\mathrm{Tr}_{\mathrm{B}}\left\{ \tilde{I}\mathcal{L}\rho_{\mathrm{B}}\right\} \left|m\right\rangle \left\langle l\right|\right\} \end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray} &&\fl \hphantom{\left[\mathcal{L}_{\iota}\right]_{ij,ml} } = \sum_{q\in L}\left|t_{q}\right|^{2}\mathrm{Tr}\left\{ \left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}dd^{\dagger}a_{q}^{\dagger}a_{q}\rho_{\mathrm{B}}\left|m\right\rangle \left\langle l\right|\right\} -\sum_{q\in L}\left|t_{q}\right|^{2}\mathrm{Tr}\left\{ \left(\left|i\right\rangle \left\langle j\right|\right)^{\dagger}da_{q}^{\dagger}\rho_{\mathrm{B}}\left|m\right\rangle \left\langle l\right|a_{q}d^{\dagger}\right\}\!\! \nonumber\\ \end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray} &&\fl \hphantom{\left[\mathcal{L}_{\iota}\right]_{ij,ml} } = \sum_{q\in L}\left|t_{q}\right|^{2}\mathrm{Tr}_{\mathrm{B}}\left\{ \left\langle i\right|dd^{\dagger}f_{q}\rho_{\mathrm{B}}\left|m\right\rangle \left.\left\langle l\right|j\right\rangle \right\} -\sum_{q\in L}\left|t_{q}\right|^{2}\mathrm{Tr}_{\mathrm{B}}\left\{ \left\langle i\right|da_{q}^{\dagger}\rho_{\mathrm{B}}\left|m\right\rangle \left\langle l\right|a_{q}d^{\dagger}\left|j\right\rangle \right\} , \nonumber\\ \end{eqnarray} \tag{ 81 }$$

where $f_{q}=\frac {1}{1+\mathrm {e}^{\beta (\varepsilon _{q}-\mu _{L})}}$ is the Fermi–Dirac distribution. In the above equation, we have used the factorized initial conditions by assuming an equilibrium Fermi distribution at t = 0 in the baths (it is worth noting that this distribution is allowed to evolve freely under the influence of the full Hamiltonian in the reduced dynamics formalism, yet the full details of bath dynamics are no longer accessible from the information stored in σ(t)). With this, it is straightforward to show that

$$\begin{eqnarray} \left[\mathcal{L}_{\iota}\right]_{ij,ml} & = & \delta_{i0}\delta_{m0}\delta_{lj}\Delta^{<}(0)-\delta_{i0}\delta_{j0}\delta_{m1}\delta_{l1}\Delta^{>}(0), \end{eqnarray} \tag{ 82 }$$

where the hybridization functions

$$\begin{eqnarray} \Delta^{<}(t) & = & \int\mathrm{d}\omega\,{\rm e}^{{\rm i}\omega t}\sum_{q\in L}\left|t_{q}\right|^{2}f_{q}\delta\left(\omega-\omega_{q}\right), \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} \Delta^{>}(t) & = & \int\mathrm{d}\omega\,{\rm e}^{{\rm i}\omega t}\sum_{q\in L}\left|t_{q}\right|^{2}\left(1-f_{q}\right)\delta\left(\omega-\omega_{q}\right) \end{eqnarray} \tag{ 84 }$$

can easily be evaluated in terms of the frequency–space coupling density given in equation (50).

The final object we need to evaluate is Φ_ι. Since no new conceptual issues arise here as compared with the calculation performed for Φ, we will simply write down the final answer. With the definitions

$$\begin{eqnarray} &&\fl \varphi_{\iota,1}^{<} \equiv \sum_{q\in L,q^{\prime}}t_{q}t_{q^{\prime}}^{*}\left\langle a_{q}^{\dagger}(t)a_{q^{\prime}}(t)\right\rangle -\sum_{q\in L}\left(\left|t_{q}\right|^{2}\left\langle d^{\dagger}(t)d(t)\right\rangle +t_{q}\varepsilon_{q}\left\langle d(t)a_{q}^{\dagger}(t)\right\rangle \right), \end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray} &&\fl \varphi_{\iota,2} \equiv \sum_{q\in L}t_{q}\left\langle d(t)a_{q}^{\dagger}(t)\right\rangle , \end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray} &&\fl \psi_{\iota,1} \equiv \sum_{q\in L,q^{\prime}}t_{q}t_{q^{\prime}}^{*}\left\langle d^{\dagger}(t)a_{q}^{\dagger}(t)a_{q^{\prime}}(t)\right\rangle , \end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray} &&\fl \psi_{\iota,2} \equiv \sum_{q\in L,q^{\prime}}t_{q}t_{q^{\prime}}\left\langle da_{q^{\prime}}^{\dagger}(t)a_{q}^{\dagger}(t)\right\rangle , \end{eqnarray} \tag{ 88 }$$

Φ_ι(t) takes the simple form

$$\begin{eqnarray} \left[\Phi_{\iota}\right]_{ij,kl}(t) & = & \left\langle l\right|\left\{ \delta_{i0}\delta_{j0}\left(\varepsilon\varphi_{\iota,2}+\varphi_{\iota,1}^{<}\right) +\delta_{i0}\delta_{j1}\left(\psi_{\iota,1}-\psi_{\iota,2}\right)\right\} \left|k\right\rangle . \end{eqnarray} \tag{ 89 }$$

The inherent asymmetry of the expression above is due to the asymmetric definition of $\tilde {I}$ (a symmetric definition would have generated non-zero matrix elements at i = j = 1 and at i = 1, j = 0, and thus our choice was motivated by computational economy). As before, it is easy to show that the block structure is such that the current can be determined without reference to the off-diagonal elements of either σ(t) or ι(t). Furthermore, it is straightforward to show that for the resonant level model (and for the Holstein model) $\left [\Phi _{\iota }\right ]_{00,00}(t)=\frac {\mathrm {d}\tilde {I}}{\mathrm {d}t}(0)$ for an initially empty dot and $\left [\Phi _{\iota }\right ]_{00,11}(t)=\frac {\mathrm {d}\tilde {I}}{\mathrm {d}t}(1)$ for an initially occupied dot. These are useful relations as they provide an alternative way of computing κ_ι(t) directly from the left current (at short times), without the need to evaluate φ_ι,2(t) or φ^<_ι,1(t).

The physics of the resonant level model are generally well known, yet the literature has seen little exploration of the properties of the memory kernel in this model, and of course none of the current memory kernel which has been introduced here. We therefore present some results below which we expect to be of interest to the field, as they provide insight into those aspects of the problem which do not rely on interaction. In order to restrict the parameter space explored, we will discuss the symmetric case in which ε = 0 with a bias voltage applied symmetrically such that V =2μ_L = −2μ_R (from here on we set ℏ = e = 1). The lead coupling densities are taken to be $\Gamma _{\mathrm {L,R}}\left (\omega \right )=\frac {1}{ (1+\mathrm {e}^{\nu (\omega -\Omega _{\mathrm { C}})} )(1+\mathrm {e}^{\nu (-\omega -\Omega _{\mathrm { C}})})}$. We will limit our attention only to the diagonal elements of the reduced density matrix and the corresponding element of the memory kernel; these elements are completely decoupled from the off-diagonal coherences, as discussed above, and therefore no approximation ensues from this.

All the results presented in this section are exact and have been calculated by a direct solution of the full equation of motion of the complete density matrix, a technique which relies on the quadratic form of the Hamiltonian and is therefore applicable only to the non-interacting case. In general, making similar progress for interacting systems requires a numerical solver of one type or another [25, 27, 29, 30, 34].

We begin with a discussion of the behavior of the memory kernel. The symmetrical parameters we have chosen are of particular interest because in the absence of interaction both σ(t) and κ(t) are completely independent of both the temperature and voltage. In addition, all non-zero matrix elements of κ are all identical to within a sign (κ_00,00 = κ_11,11 = −κ_11,00 = −κ_00,11). In figure 1 we therefore explore the dependence of one arbitrarily chosen element of κ on a range of band parameters: cutoff energies Ω_C (the bandwidth is ∼2Ω_C) and band cutoff widths $\frac {1}{\nu }$. In each panel we go from a small bandwidth (red) to a large one (blue) at a set cutoff width, with the sharpest cutoff shown in panel A, an intermediate value in B and the smoothest in C.

Zoom In Zoom Out Reset image size

The effect of the two parameters describing our chosen band shape on the memory kernel can be understood quite well by considering the trends shown in the plot: as ν decreases, reflections are softened by the gradual slope at the band edge and the memory kernel decays more quickly. On the other hand, increasing the bandwidth induces oscillations at a frequency ω ≈ Ω_C, but also increases the proportional weight of the short-time part of the memory kernel. Eventually, if we were to approach the wide band limit, the memory would approach the form of a delta function and a Markovian description of the dynamics would become exact.

The current memory kernel κ_ι is not as highly symmetric as κ, and depends to some extent on all the parameters of the problem. There are two distinct (though similar) elements in κ_ι at non-zero voltage, and these are plotted for the left current with Γν = 10 in panels A and B of figure 2. The figure illustrates the dependence of κ_ι on the bandwidth, which exhibits the same properties observed in κ. This is also true of its ν dependence: to exemplify this, figure 3 displays the same data for Γν = 0.5, where κ_ι decays more quickly and smoothly. Interestingly, the timescale over which κ_ι decays to zero does not appear to differ markedly from the corresponding timescale for κ at similar parameters; this suggests that the cutoff approximation remains as useful for the current as it is for impurity observables.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effect of temperature on κ_ι depends greatly on the choice of other parameters. At the parameters we have chosen for figure 4 the asymmetry between the two κ_ι elements is increased somewhat when the temperature is lowered, corresponding to an increase in the current (not shown). One would expect a rather different effect when, for instance, things are set up in such a way that thermal enhancement of the current occurs. Unlike κ, κ_ι depends on temperature even in the fully symmetric and non-interacting case is interesting, and expresses the fact that this quantity is connected to bath observables as well as to those in the impurity subspace.

Zoom In Zoom Out Reset image size

Figure 5 shows how voltage affects the current memory kernel: at zero voltage the two elements of κ_ι are identical up to a sign, and the application of a voltage increases the diagonal element while suppressing the off-diagonal terms. Additionally, an increase in the oscillation frequency is observed for κ_00,11, but no such clear trend exists for the oscillation frequency in κ_00,00. As V passes the bandwidth (∼20Γ here), the left lead becomes entirely occupied and the right entirely empty, and further increasing the voltage ceases to have any effect on the current memory kernel, just as occurs in the case of the current itself (not shown).

Zoom In Zoom Out Reset image size

To show that the generalized NZME formalism introduced in this work indeed reproduces the correct results for the current as a function of time, figure 6 presents a comparison between the current obtained directly (lighter solid lines) and by way of equation (29) (darker dashed lines). Pairs of lines describing the two different ways of obtaining currents at identical parameters overlap to within numerical errors, expressing the equivalence between the two methods when convergence in the cutoff time has been attained and the correctness of the approach at the t_c → ∞ limit.

Zoom In Zoom Out Reset image size

Finally, while the NZME memory kernel technique and the cutoff approximation have been shown to be efficient for σ in a variety of interacting and non-interacting cases [4, 36, 40], meaning that results at long times t ≫ t_c converge at a finite t_c, no such calculations have previously been carried out for the generalized technique. This entails a convergence analysis of the type exemplified graphically in figure 7. In essence, the cutoff time t_c must be increased until the desired accuracy is reached. In the example shown in figure 7, convergence is achieved quickly and even short-time oscillations beyond the range of t_c are predicted with some accuracy (as can be seen from the extension of the oscillatory ridges beyond the boundary of the transparent t = t_c plane). As an alternative (if partial) representation of this idea, in figure 8 several plane cuts through this function are shown, but this time at Γν = 0.5; both the current and its time derivative are displayed. The rapid convergence visible in either representation illustrates that the idea of the reduced dynamics technique and the cutoff approximation remains useful in practice even for the current, despite it being a non-system operator not accessible within the confines of the standard NZME formalism.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have reviewed the process of implementing reduced dynamics techniques by way of the NZME and the memory cutoff approximation, which have recently been introduced with great success into several numerically exact non-equilibrium impurity solvers. The procedure of deriving a memory kernel scheme for a general impurity model and obtaining calculable expressions in terms of standard second-quantization operators was outlined, and the example of the non-interacting resonant level model was worked out in full detail. For this non-interacting case, some illustrative examples of the physical properties of the memory kernel were discussed.

An important limitation of the reduced dynamics techniques so far has been the lack of access to non-impurity observables, such as the electronic current in the resonant level model, the Anderson impurity model and the Holstein model. A generalization of the NZME formalism which allows access to general operators was therefore introduced here, and the implementation of this formalism was carried through for the example of the current in the non-interacting limit. This led to the definition and evaluation of a current memory kernel κ_ι, which was subsequently explored for its dependence on time, bandwidth, voltage and temperature. The validity of the cutoff approximation for the current memory kernel was then verified and discussed.

Looking forward, we expect the ideas expounded upon here to have several major implications: firstly, we hope to see them become a standard part of the toolbox of high quality time domain numerical simulations of impurity models, and extended to a variety of models and methods; in this context the memory technique should be viewed not as competing with existing direct solvers, but as a supplemental tool which allows efficient extension of any general short-time solver to long timescales, in situations where the memory timescale is short. Secondly, since access to current enables access to Green's functions, the benefits offered by memory techniques are expected to be applicable to interacting lattice simulations as well, by way of mapping schemes such as dynamical mean field theory and its various extensions. Finally, we believe the memory kernel framework is a fertile ground for defining new approximation schemes more general than the cutoff approximation, and in this context it will be particularly interesting to understand the long-time behavior of the memory kernel in interacting cases and its behavior in larger impurity models.

We thank A Nitzan and M R Wegewijs for insightful comments and helpful conversations. GC is grateful to Yad Hanadiv–Rothschild Foundation for the award of a Rothschild Postdoctoral Fellowship. EYW is grateful to The Center for Nanoscience and Nanotechnology at Tel Aviv University for a doctoral fellowship. This work was supported by the US–Israel Binational Science Foundation.