Lindblad equation and its semi-classical limit of the Anderson-Holstein model

TL;DR

This paper investigates the approximation of the Anderson-Holstein model, a multi-level open quantum system, using Redfield and Lindblad equations, and derives their semi-classical limits with explicit classical master equations.

Contribution

It provides a detailed comparison of Redfield and Lindblad equations for the Anderson-Holstein model and derives their semi-classical limits explicitly.

Findings

Semi-classical limits of Redfield and Lindblad equations are explicitly derived.

Classical master equations share similar transport and hopping structures.

Comparison of Redfield and Lindblad equations via time-dependent perturbation theory.

Abstract

For multi-level open quantum system, the interaction between different levels could pose challenge to understand the quantum system both analytically and numerically. In this work, we study the approximation of the dynamics of the Anderson-Holstein model, as a model of multi-level open quantum system, by Redfield and Lindblad equations. Both equations have a desirable property that if the density operators for different levels is diagonal initially, they remain to be diagonal for any time. Thanks to this nice property, the semi-classical limit of both Redfield and Lindblad equations could be derived explicitly; the resulting classical master equations share similar structures of transport and hopping terms. The Redfield and Lindblad equations are also compared from the angle of time dependent perturbation theory.

Equations

$$\begin{array}{rl}\hat{H}& ={\hat{H}}_s+{\hat{H}}_b+{\hat{H}}_c\\ {\hat{H}}_s& =\frac{{\hat{p}}^2}{2m}+\frac{1}{2}m{\omega }_s^2{\hat{x}}^2+U(\hat{x}){\hat{d}}^{†}\hat{d}\\ {\hat{H}}_b& =\sum _k(E_k-\mu ){\hat{c}}_k^{†}{\hat{c}}_k\\ {\hat{H}}_c& =\sum _k{V}_k({\hat{c}}_k^{†}\hat{d}+{\hat{c}}_k{\hat{d}}^{†})\end{array}$$

$$ϵ:=\frac{\hslash {\omega }_s}{\mathsf{E}}$$

$$\alpha :=\frac{\mathsf{V}}{\mathsf{E}}$$

$$\frac{V_k}{\mathsf{E}}=\alpha {\tilde{V}}_k.$$

$${\hat{H}}_{non}$$

$${\hat{H}}_{s,non}$$

$${\hat{H}}_{b,non}$$

$${\hat{H}}_{c,non}$$

$$i\epsilon\partial_{\tilde{t}}\hat{\rho}=\mathopen{\big{missing}}[\hat{H}_{non},\hat{\rho}\mathopen{\big{missing}}],$$

$$\begin{array}{rl}\hat{H}& ={\hat{H}}_s+{\hat{H}}_b+\alpha {\hat{H}}_c\\ {\hat{H}}_s& =-\frac{1}{2}{ϵ}^2{\nabla }_x^2+\frac{1}{2}{x}^2+U(x){\hat{d}}^{†}\hat{d}\\ {\hat{H}}_b& =\sum _k{E}_k{\hat{c}}_k^{†}{\hat{c}}_k\\ {\hat{H}}_c& =\sum _k{V}_k({\hat{c}}_k^{†}\hat{d}+{\hat{c}}_k{\hat{d}}^{†})=:{\hat{C}}^{†}\hat{d}+\hat{C}{\hat{d}}^{†}\end{array}$$

$$U(x):=\sqrt{2}gx+g^2+{\bar{U}}_0,$$

$$\hat{C}:=\sum _k{V}_k{\hat{c}}_k.$$

$$i\epsilon\partial_{t}\hat{\rho}=\mathopen{\big{missing}}[\hat{H},\hat{\rho}\mathopen{\big{missing}}].$$

$${\hat{H}}_0:=⟨0|{\hat{H}}_s|0⟩\equiv -\frac{1}{2}{ϵ}^2{\nabla }_x^2+\frac{1}{2}{x}^2$$

$${\varphi }_k^0=N_k^0{H}_k\left(\frac{x}{\sqrt{ϵ}}\right)\exp \left(-\frac{x^2}{2ϵ}\right),k\in \mathbb{N},$$

$${\hat{H}}_1:=⟨1|{\hat{H}}_s|1⟩\equiv -\frac{1}{2}{ϵ}^2{\nabla }_x^2+\frac{1}{2}(x+\sqrt{2}g{)}^2+{\bar{U}}_0$$

$${\hat{O}}_I(t):=e^{\frac{i}{ϵ}({\hat{H}}_s+{\hat{H}}_b)t}\hat{O}{e}^{-\frac{i}{ϵ}({\hat{H}}_s+{\hat{H}}_b)t}.$$

$$\frac{d}{dt}\hat{\rho}_{I}(t)=-i\frac{\alpha}{\epsilon}\mathopen{\big{missing}}[\hat{H}_{c,I}(t),\hat{\rho}_{I}(t)\mathopen{\big{missing}}]=:\eta\mathcal{L}(t)\hat{\rho}_{I}(t)$$

$${\hat{H}}_{c,I}(t)={\hat{C}}_I^{†}(t){\hat{d}}_I(t)+{\hat{d}}_I^{†}(t){\hat{C}}_I(t),$$

$${\hat{d}}_I(t)=e^{\frac{i}{ϵ}{\hat{H}}_{s}t}\hat{d}{e}^{-\frac{i}{ϵ}{\hat{H}}_{s}t},\text{and}{\hat{C}}_I(t)=\sum _k{V}_k{e}^{-\frac{i}{ϵ}{E}_{k}t}{\hat{c}}_k.$$

$$\mathcal{P}A={\text{Tr}}_b(A)\otimes {\hat{\rho }}_{b,eq}$$

$${\hat{\rho }}_I(s)=G(t,s){\hat{\rho }}_I(t)=G(t,s)(\mathcal{P}+\mathcal{Q}){\hat{\rho }}_I(t),(\text{for }s\le t)$$

$$\frac{d}{dt}\mathcal{Q}{\hat{\rho }}_I(t)=\eta \mathcal{Q}\mathcal{L}(t)\mathcal{P}{\hat{\rho }}_I(t)+\eta \mathcal{Q}\mathcal{L}(t)\mathcal{Q}{\hat{\rho }}_I(t).$$

$$\begin{array}{rl}\mathcal{Q}{\hat{\rho }}_I(t)& =\mathcal{G}(t,t_0)\mathcal{Q}{\hat{\rho }}_I(t_0)+\eta {\int }_{t_0}^{t}ds\mathcal{G}(t,s)\mathcal{Q}\mathcal{L}(s)\mathcal{P}{\hat{\rho }}_I(s)\\ & \stackrel{\text{eqref}eq:solvonneu}{=}\mathcal{G}(t,t_0)\mathcal{Q}{\hat{\rho }}_I(t_0)+\eta {\int }_{t_0}^{t}ds\mathcal{G}(t,s)\mathcal{Q}\mathcal{L}(s)\mathcal{P}G(t,s)(\mathcal{P}+\mathcal{Q}){\hat{\rho }}_I(t)\\ & =\mathcal{G}(t,t_0)\mathcal{Q}{\hat{\rho }}_I(t_0)+\Sigma (t)(\mathcal{P}+\mathcal{Q}){\hat{\rho }}_I(t)\end{array}$$

$$\mathcal{Q}{\hat{\rho }}_I(t)={\left(\text{Id}-\Sigma (t)\right)}^{-1}\Sigma (t)\mathcal{P}{\hat{\rho }}_I(t)$$

$$\begin{array}{rl}\frac{d}{dt}\mathcal{P}{\hat{\rho }}_I(t)& =\eta \mathcal{P}\mathcal{L}(t)\mathcal{P}{\hat{\rho }}_I(t)+\eta \mathcal{P}\mathcal{L}(t)\mathcal{Q}{\hat{\rho }}_I(t)\\ & =\eta \mathcal{P}\mathcal{L}(t)\mathcal{P}{\hat{\rho }}_I(t)+\eta \mathcal{P}\mathcal{L}(t)(\text{Id}-\Sigma (t){)}^{-1}\Sigma (t)\mathcal{P}{\hat{\rho }}_I(t)\\ & =\eta \mathcal{P}\mathcal{L}(t)(\text{Id}-\Sigma (t){)}^{-1}\mathcal{P}{\hat{\rho }}_I(t)\\ & =\left(\eta \mathcal{P}\mathcal{L}(t)\mathcal{P}+{\eta }^2\mathcal{P}\mathcal{L}(t){\int }_{t_0}^{t}ds\mathcal{Q}\mathcal{L}(s)\mathcal{P}\right){\hat{\rho }}_I(t)+\mathcal{O}({\eta }^3)\end{array}$$

$$\begin{array}{rl}\frac{d}{dt}\mathcal{P}{\hat{\rho }}_I(t)& ={\eta }^2{\int }_{t_0}^{t}ds\mathcal{P}\mathcal{L}(t)\mathcal{L}(s)\mathcal{P}{\hat{\rho }}_I(t)\end{array}$$

$$\frac{d}{dt}\hat{\rho}_{s,I}(t)=-\left(\frac{\alpha}{\epsilon}\right)^{2}\text{Tr}_{b}\left(\int_{t_{0}}^{t}\ ds\mathopen{\big{missing}}[\hat{H}_{c,I}(t),[\hat{H}_{c,I}(s),\hat{\rho}_{s,I}(t)\otimes\hat{\rho}_{b,eq}]\mathopen{\big{missing}}]\right).$$

$$\frac{d}{dt}{\hat{\rho }}_{s,I}(t)=-{\left(\frac{\alpha }{ϵ}\right)}^2{\text{Tr}}_b\left({\int }_0^{\infty }d\tau \left[{\hat{H}}_{c,I}(t),[{\hat{H}}_{c,I}(t-\tau ),{\hat{\rho }}_{s,I}(t)\otimes {\hat{\rho }}_{b,eq}]\right]\right)$$

$$\begin{split}\frac{d}{dt}\hat{\rho}_{s,I}(t)=-\left(\frac{\alpha}{\epsilon}\right)^{2}\int_{0}^{\infty}d\tau\ &\mathopen{\big{missing}}[\hat{d}_{I}(t),\hat{d}^{\dagger}_{I}(t-\tau)\hat{\rho}_{s,I}(t)\mathopen{\big{missing}}]F(t,t-\tau)\\ +&\mathopen{\big{missing}}[\hat{\rho}_{s,I}(t)\hat{d}^{\dagger}_{I}(t-\tau),\hat{d}_{I}(t)\mathopen{\big{missing}}]G(t,t-\tau)+h.c.\\ \end{split}$$

$$\begin{array}{rl}F(t,t^{\prime })& :={\text{Tr}}_b\left({\hat{C}}_I^{†}(t){\hat{C}}_I(t^{\prime }){\hat{\rho }}_{b,eq}\right)=\sum _k{V}_k^2\exp \left(\frac{i}{ϵ}{E}_k(t-t^{\prime })\right)f(E_k)\\ G(t,t^{\prime })& :={\text{Tr}}_b\left({\hat{C}}_I(t^{\prime }){\hat{C}}_I^{†}(t){\hat{\rho }}_{b,eq}\right)=\sum _k{V}_k^2\exp \left(\frac{i}{ϵ}{E}_k(t-t^{\prime })\right)(1-f(E_k))\end{array}$$

$$\begin{split}\frac{d}{dt}\hat{\rho}_{s}(t)=-\frac{i}{\epsilon}\mathopen{\big{missing}}[\hat{H}_{s},\hat{\rho}_{s}(t)\mathopen{\big{missing}}]-\left(\frac{\alpha}{\epsilon}\right)^{2}\int_{0}^{\infty}d\tau\ &\left(\hat{d}e^{-\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{d}^{\dagger}e^{\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{\rho}_{s}(t)-e^{-\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{d}^{\dagger}e^{\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{\rho}_{s}(t)\hat{d}\right)F(t,t-\tau)\\ +&\left(\hat{\rho}_{s}(t)e^{-\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{d}^{\dagger}e^{\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{d}-\hat{d}\hat{\rho}_{s}(t)e^{-\frac{i}{\epsilon}\hat{H}_{s}\tau}\hat{d}^{\dagger}e^{\frac{i}{\epsilon}\hat{H}_{s}\tau}\right)G(t,t-\tau)\\ +&h.c.\\ \end{split}$$

Full text

Multi-level open quantum system, Anderson-Holstein model, Redfield equation, Lindblad equation, semi-classical limit

I Introduction

Multi-level open quantum systems have received much attention due to their wide applications and intriguing phenomena Ghosh et al. (2004); Galperin et al. (2007a, 2008). One of the simplest models is perhaps the Anderson-Holstein model, a two-level open quantum system Holstein (1959). The Anderson-Holstein model is a simplistic model for a molecule as the system of interest, represented by a classical nucleus degree of freedom and a two level electronic degree of freedom, coupled with a bath of fermions, for instance, a reservoir of electrons. In this paper, the simplified version of Anderson-Holstein model discussed in Dou et al. (2015a) is used as an example for illustrating purpose, which will be explained in more details in the next section. The goal is to understand the approach of quantum master equations, in particular, the Lindbladian formalism, for such systems in the weak-coupling limit and also to study the semiclassical limit of the quantum master equations. Our study here should generalize to other multi-level open quantum systems.

Anderson-Holstein model, since introduced, has been widely studied using various theoretical and numerical approaches, for instance, the Green’s function approach Ryndyk (2016); Galperin et al. (2004), equation-of-motion method Galperin et al. (2006, 2007b), quantum Monte Carlo method Mühlbacher and Rabani (2008), semi-classical approximation Mitra et al. (2004, 2005), non-crossing approximation Chen et al. (2016) and by using quantum master equations Elste et al. (2008); Esposito and Galperin (2009, 2010); Dou et al. (2015a, b, c); Dou and Subotnik (2016). In the perspective of quantum master equation, which is mostly related to the current work, the quantum master equation in Redfield flavor for Anderson-Holstein model has been derived in Elste et al. (2008); Esposito and Galperin (2009, 2010). The semiclassical limit of the Redfield equation, known as the classical master equation (CME), has been considered in Dou et al. (2015a), which also proposed a numerical method based on surface hopping. The CME perspective has then been used to study various physical aspects of Anderson-Holstein model, e.g., broadening, Marcus rate Dou et al. (2015b, c); Dou and Subotnik (2016). In all these works, the focus has been on Redfield equation (or Redfield generator). As far as we know, not much attention has been put into the Lindbladian formulation nor its semiclassical limit of the Anderson-Holstein model, which is the focus of the current work.

Recall that closed quantum systems can be fully characterized by the Hamiltonian; its time-evolution dynamics is characterized by Schrödinger equation (or von Neumann equation if we are dealing with density operators). In the framework of quantum master equation, open quantum systems can be described by Nakajima-Zwanzig equation with the assistance of projection operator to a subspace in which density operator for closed system is separable Breuer and Petruccione (2002). Although Nakajima-Zwanzig equation provides an exact expression for the open quantum system in the interaction picture, in general, it is not easy to retrieve useful information, neither analytically nor numerically. Part of the reason is attributed to the memory effect in the Nakajima-Zwanzig equation. While it is, of course, important to research on non-Markovian dynamics itself; many questions in non-Markovian dynamics are still open and the mathematical foundation requires further investigation Breuer et al. (2016). Often times Markovian approximation is taken to simplify the governing equations, which is a valid approximation in the weak coupling regime. See also Davies (1974) for mathematical study of the quantum Markovian approximation.

The Markovian approximation leads to Redfield equation using time-convolutionless (TCL) projection operator method in the weak-coupling limit Breuer and Petruccione (2002). Furthermore, with secular approximation, Lindblad equation can be obtained from Redfield equation Breuer and Petruccione (2002); the Lindblad equation has better mathematical structures such as complete positivity Lindblad (1976). The essence of secular approximation is to approximate fast oscillating terms by zero in the sense of averaging on a coarser time scale. The goal of this paper is to study Redfield equation and Lindblad equation for multi-level open quantum systems as well as their semi-classical limit, in the context of the Anderson-Holstein model. For these two equations, one obtains the classical master equations (CME) and Lindbladian classical master equations (LCME) in the semiclassical limit. The relations of various models and the asymptotic limit connecting those are summarized in Figure 1.

It is worth mentioning that there is a debate in the literature on which equation better models the open quantum system, especially when the coupling between the system and the bath is not weak. The underlying discussion focuses on whether complete positivity (CP) is necessary for modeling open quantum systems. There are at least two arguments supporting complete positivity in quantum systems: the first one is from the perspective of “total domain”; the second one from “product state” McCracken (2013). What’s more, one recent research indicates that without complete positivity in Redfield equation, the dynamics might be inconsistent with the second law of thermodynamics Argentieri et al. (2014). Some, however, criticize that we might over-emphasize the importance of complete positivity in modeling open quantum systems. Pechukas proposed that for a composite quantum system with entangled initial condition, the positivity property might not hold for the reduced dynamics Pechukas (1994). Shaji and Sudarshan argued that complete positivity is not necessary by carefully examinizing arguments supporting complete positivity Shaji and Sudarshan (2005). Negativity, as opposed to positivity, is not only observed in experiment but also can be informative to the coupling with bath McCracken (2013). In our study on Anderson-Holstein model, imposing complete positivity (and thus Lindblad equation) should be justified as we only consider the weak-coupling regime. In particular, as will become clear in our analysis, under the same assumption used in deriving Redfield equation, the secular approximation for getting Lindblad equation is in fact also justified; hence, the use of Lindblad equation is natural.

In this paper, we consider Anderson-Holstein model in weak-coupling and semi-classical limits. Under the assumption that the coupling strength is weak, we will revisit the derivation of Redfield equation in Section III and derive the explicit form of Lindblad equation in Section IV. The semi-classical study of both equations is discussed in Section V. The perturbation result of both equations is presented in Section VI. Section VII summarizes the main results and ends the paper with some concluding remarks.

II Anderson-Holstein model

The Anderson-Holstein model under study here describes a two-level system coupled with a bath of many non-interacting electrons (or in general spin-$1/2$ fermions). For instance, the two-level system can be thought as a simplistic model for the nuclei degree of freedom of a molecule with two potential energy surfaces depending on the electronic state of the molecule. For simplicity, in the Anderson-Holstein model, the two-level system is in one spatial dimension and one of the potential energy surface is taken to be a harmonic oscillator with frequency $\omega_{s}$, and the difference $U(x)$ between the two potential energy surface is modeled as a linear function of the nucleus position (and thus is also a harmonic oscillator with shifted center and energy) Holstein (1959). More specifically, the Hamiltonian for the whole system is given by

[TABLE]

where we follow the notation of Dou et al. (2015a): $\hat{d}$ and $\hat{d}^{\dagger}$ are the annihilation and creation operators for the two-level electron state of the molecule, $\hat{c}_{k}$ and $\hat{c}^{\dagger}_{k}$ are the annihilation and creation operators for electron states in the bath, $E_{k}$ is the energy level of those states, $\mu$ is the Fermi level, and $V_{k}$ is the coupling strength between the molecule and the $k$-th mode in bath, assumed to be real. The Hilbert space corresponds to the molecule is thus L^{2}(\mathbb{R})\otimes\mathbb{C}^{2}=L^{2}(\mathbb{R})\otimes\operatorname{span}\mathopen{\big{missing}}\{\lvert 0\rangle,\lvert 1\rangle\mathclose{\big{missing}}\}, so that we have $\hat{d}\lvert 1\rangle=\lvert 0\rangle$ and $\hat{d}^{\dagger}\lvert 0\rangle=\lvert 1\rangle$. Thus in (II.1), $\hat{H}_{s}$ is the Hamiltonian operator of the “system”, $\hat{H}_{b}$ is the Hamiltonian of the “bath”, and $\hat{H}_{c}$ describes the coupling between the system and the bath. The goal is to understand the evolution of the system as an open quantum system (i.e., integrate out the bath degree of freedom). The above is a simplified version of the original Anderson-Holstein model as (1) only one electron is considered for the molecule, thus the model excludes Coulomb interaction between electrons of the molecule; (2) the molecule is coupled to one electrode, but not two electrodes, and thus only one Fermi level $\mu$ is used for the environment.

Remark*.*

In the literature, e.g. Mitra et al. (2005), sometimes the term “single-level” is used for the above system to emphasize that there is only one on-site electron. We have adopted here the term “two-level” to emphasize that the Hilbert space for molecular system is L^{2}(\mathbb{R})\otimes\text{span}\mathopen{\big{missing}}\{\lvert 0\rangle,\lvert 1\rangle\mathclose{\big{missing}}\} and the second component has dimension $2$.

To proceed, let us first non-dimensionalize the problem according to the following rescaling:

Denote $\mathsf{\ell}$ the characteristic length scale of $x$, the position degree of freedom of the system. That is, if we take $x=\mathsf{\ell}\tilde{x}$, $\tilde{x}$ becomes a dimensionless quantity with order $\mathcal{O}(1)$. 2. 2.

As a consequence, the scaling factor for molecular energy is then $\mathsf{E}=m\omega_{s}^{2}\mathsf{\ell}^{2}$. We will use $\mathsf{E}$ as the scaling factor for all physical quantities whose dimension is energy (thus including all terms in the Hamiltonian). 3. 3.

Denote $\mathsf{T}$ the time scale of the evolution of the system. Thus, $t=\mathsf{T}\tilde{t}$ where $\tilde{t}=\mathcal{O}(1)$. Physically, it is reasonable to choose $\mathsf{T}=\frac{1}{\omega_{s}}$ Shankar (1994), since this is the time scale of an isolated harmonic oscillator with frequency $\omega_{s}$. 4. 4.

As both $\hbar\omega_{s}$ and $\mathsf{E}$ have the energy dimension, the ratio

[TABLE]

is a dimensionless quantity. In our analysis of the semi-classical limit of the system, we will assume that $\epsilon$ is a small parameter $\epsilon\downarrow 0$. 5. 5.

Let $\mathsf{V}$ be the typical interaction strength with dimension as energy, i.e., we assume $\tilde{V}_{k}:=\frac{V_{k}}{\mathsf{V}}=\mathcal{O}(1)$. The ratio of $\mathsf{V}$ and $\mathsf{E}$

[TABLE]

is dimensionless. Weak-coupling limit means $\mathsf{V}\ll\mathsf{E}$, and hence $\alpha\downarrow 0$, which means physically that the coupling between the bath and the system is weak compared to the typical energy scale of the system. Note that by the above rescaling, we have

[TABLE]

After performing the above rescaling, the non-dimensionalized Hamiltonian becomes

[TABLE]

where $\tilde{U}(\tilde{x})=\frac{U(x)}{\mathsf{E}}$, $\tilde{E}_{k}=\frac{E_{k}}{\mathsf{E}}$ and $\tilde{\mu}=\frac{\mu}{\mathsf{E}}$. Moreover, the von Neumann equation becomes

[TABLE]

where $\hat{\rho}$ is the density operator for the closed system and \mathopen{\big{missing}}[\ ,\ \mathopen{\big{missing}}] is the usual commutator.

Dropping “tilde” and “non” and shifting the energy reference to replace $E_{k}-\mu$ by $E_{k}$ to simplify the notation, we arrive at the Hamiltonian

[TABLE]

Here, $U(x)$ is a linear function of position $x$ with the form (the reason of the specific choice or parametrization will become clear below)

[TABLE]

where $\bar{U}_{0}$ is known as renormalized energy, and the weighted annihilation operator $\hat{C}$ is defined as

[TABLE]

The evolution of the density operator is given by the von Neumann equation

[TABLE]

All quantities in the Equations (II.2) and (II.3) are dimensionless and the parameters $V_{k}$ and $E_{k}$’s are $\mathcal{O}(1)$.

In summary, after non-dimensionalization, the model contains two scaling parameters $\epsilon$ and $\alpha$, corresponding to the semiclassical parameter and coupling strength respectively. In the rest of the paper, we will consider the weak-coupling limit and the semiclassical limit. In the weak-coupling limit, we have $\alpha\downarrow 0$, which leads to Redfield and Lindblad equations for fixed $\epsilon$, while the semiclassical limit means $\epsilon\downarrow 0$. See Figure 1 for an overview.

Notice that

[TABLE]

is a Hamiltonian for a single harmonic oscillator. It is well-known that it has eigenfunctions

[TABLE]

where $H_{k}$ is the $k$-th Hermite polynomial and $N_{k}^{0}$ is a normalization constant. The corresponding eigenvalue is $\epsilon(k+\frac{1}{2})$. Similarly,

[TABLE]

has eigenfunctions $\phi^{1}_{k}=\phi_{k}^{0}(x+\sqrt{2}g)$ and corresponding eigenvalues $\epsilon(k+\frac{1}{2})+\bar{U}_{0}$. Note that this explains the specific parametrization of the linear function $U(x)$ above. To keep notations and calculations simple, we will only consider the Anderson-Holstein model with renormalized energy $\bar{U}_{0}=0$. The extension to the general case is straightforward.

III Revisiting the Derivation of Redfield Equation

The derivation of Redfield equation has been well studied and presented in e.g., Breuer and Petruccione (2002); Blum (1996); Elste et al. (2008). Physically, Born-Markov approximation is the key to reduce the dynamics of the system to a Markovian dynamics. In this section, we will revisit the derivation for the Anderson-Holstein model to set the grounds of our discussion below, using the time-convolutionless equation (TCL) approach following (Breuer and Petruccione, 2002, Chapter 9); we will borrow notations from this reference as well. For fixed $\epsilon$, we consider the weak-coupling limit below, that is, $\alpha\downarrow 0$ while $\epsilon$ stays fixed, and hence $\eta:=\frac{\alpha}{\epsilon}\downarrow 0$.

To simplify the dynamical equation, it is more convenient to use the interaction picture (with respect to the uncoupled system and bath), so that the operators are given by

[TABLE]

Interaction picture is very convenient in weak-coupling limit since it removes the effect of fast motion (due to $\hat{H}_{s}+\hat{H}_{b}$) from the slow motion (due to $\alpha\hat{H}_{c}$).

The von Neumann equation in the interaction picture is

[TABLE]

where we have introduced \mathcal{L}(t):=-i\mathopen{\big{missing}}[\hat{H}_{c,I}(t),\cdot\mathopen{\big{missing}}] as a super-operator acting on density operators. By explicit calculation, we have

[TABLE]

where

[TABLE]

For any trace-class operator $A$ defined for the whole closed system, we define a projection operator $\mathcal{P}$ by

[TABLE]

where $\text{Tr}_{b}(A)$ is the partial trace over bath degree of freedom, and where $\hat{\rho}_{b,eq}:=e^{-\beta\hat{H}_{b}}/Z_{b}$ is the density operator of electron bath at thermal equilibrium, with $\beta=\frac{\mathsf{E}}{k_{B}T}$ (inverse of the rescaled temperature) and $Z_{b}$ is the partition function. This projection operator $\mathcal{P}$ disentangles the system and bath and replaces the bath by the thermal equilibrium; this is a core ingredient in Born approximation, whose physical reasoning can be found in (Blum, 1996, p 276). We also define its orthogonal complement as $\mathcal{Q}:=\text{Id}-\mathcal{P}$. For a given density matrix $\hat{\rho}_{I}(t)$, $\mathcal{P}\hat{\rho}_{I}(t)$ is known as relevant part and $\mathcal{Q}\hat{\rho}_{I}(t)$ irrelevant part.

We may formally write down the solution to Equation (III.1) using Green’s function as

[TABLE]

where $G(t,s):=T_{\rightarrow}\exp\left(-\eta\int_{s}^{t}ds^{\prime}\mathcal{L}(s^{\prime})\right)$ and $T_{\rightarrow}$ represents anti-chronological time-ordering operator. On the other hand, applying operator $\mathcal{Q}$ to Equation (III.1) gives a differential equation

[TABLE]

By Duhamel’s principle, its solution in integral form is

[TABLE]

where $t_{0}$ is the starting time of interest, $\mathcal{G}(t,s)=T_{\leftarrow}\exp\left(\eta\int_{s}^{t}ds^{\prime}\ \mathcal{Q}\mathcal{L}(s^{\prime})\right)$, $T_{\leftarrow}$ represents the chronological time-ordering operator, and $\Sigma(t):=\eta\int_{t_{0}}^{t}ds\ \mathcal{G}(t,s)\mathcal{Q}\mathcal{L}(s)\mathcal{P}G(t,s)$. Assume that at time $t_{0}$, the bath is at thermal equilibrium and density operator $\hat{\rho}_{I}(t_{0})$ is separable, i.e., $\hat{\rho}_{I}(t_{0})=\hat{\rho}_{s,I}(t_{0})\otimes\hat{\rho}_{b,eq}$. Then $\mathcal{Q}\hat{\rho}_{I}(t_{0})=0$ and hence we obtain

[TABLE]

if $\text{Id}-\Sigma(t)$ is invertible, which is the case, for instance, when $\eta$ is so small that $\lVert\Sigma(t)\rVert<1$.

We may also apply the operator $\mathcal{P}$ to Equation (III.1) and get

[TABLE]

In the last step we perform asymptotic expansion of operator in terms of $\eta$. It could be easily verified that $\mathcal{P}\mathcal{L}(t)\mathcal{P}=0$ from the definition of $\mathcal{L}$. Then the leading order expansion is

[TABLE]

After replacing $\mathcal{L}$ by its definition, we arrive at

[TABLE]

The leading order expansion of $\frac{d}{dt}\hat{\rho}_{s,I}(t)$ is the same as applying Born-Markov approximation to von Neumann equation directly. Change the variable $\tau=t-s$ and push $t_{0}$ to approach $-\infty$, then the last equation becomes Redfield equation

[TABLE]

There are two formal justifications for pushing $t_{0}$ to $-\infty$: if $t_{0}=-\infty$, the dynamics does not depend on initial time as a parameter; moreover, if the system evolves from long time ago, we may as well consider $t_{0}=-\infty$.

After opening the double commutator and simplify the equation, we arrive at

[TABLE]

where time correlation functions

[TABLE]

with $f(z)$ being the Fermi-Dirac function $f(z)=1/(1+e^{\beta z})$.

Transforming Equation (III.5) back into Schrödinger picture, we end up with

[TABLE]

This is the Redfield equation for Anderson-Holstein model.

IV Derivation of Lindblad Equation

It is a fundamental result Lindblad (1976) that a completely positive dynamical map can be written in Lindblad form for open quantum systems. Thus, in this section, we aim at deriving the Lindblad equation for Anderson-Holstein model. The derivation of Lindblad equation from a microscopic point of view has been studied for some cases, see e.g., Breuer and Petruccione (2002), though to the best of our knowledge not for the Anderson-Holstein model. In this section, we will show that under the previous condition that

[TABLE]

Lindblad equation in Schrödinger picture for the Anderson-Holstein model is given by

[TABLE]

with a Lindbladian corrected Hamiltonian

[TABLE]

and dissipative operator

[TABLE]

where $\left[\hat{A},\hat{B}\right]_{+}:=\hat{A}\hat{B}+\hat{B}\hat{A}$ is anti-commutator for two operators $\hat{A}$ and $\hat{B}$. The coefficients will be given in Equation (IV.21) below and the operators $\hat{D}^{(\dagger)}(\omega)$ will be defined below (see Equation (IV.4)). It is clear that the above dissipative operator takes the Lindblad form.

IV.1 An alternative representation for Redfield equation

In Anderson-Holstein model, it is natural to consider eigenfunctions of $\hat{H}_{s}$, which form two energy ladders. The evolution of system can be considered as quantum jumping between different energy levels. Thus the annihilation and creation operators might be decomposed in terms of the numbers of energy levels that the system jumps. Such decomposition has been used to derive Lindblad equation in (Breuer and Petruccione, 2002, p 125-131). We shall use this technique to study Redfield equation and Lindblad equation below for Anderson-Holstein model.

Definition 1**.**

For each $\omega\in\mathbb{Z}$, define an operator

[TABLE]

where $\Pi^{(m)}_{k}$ is the projection operator to quantum state $\lvert\phi^{m}_{k}\rangle\otimes\lvert m\rangle\equiv\lvert\phi_{k}^{m},m\rangle$, $m\in\{0,1\}$, $k,k^{\prime}\in\mathbb{N}$.

Recall that $\lvert\phi_{k}^{m}\rangle$ is the eigenfunction of Hamiltonian $\hat{H}_{s}$ discussed at the end of Section II. In other words, $\Pi_{k}^{(m)}:=\lvert\phi_{k}^{m},m\rangle\langle\phi_{k}^{m},m\rvert$. Hence, the adjoint operator of $\hat{D}(\omega)$ is

[TABLE]

This definition was used in Breuer and Petruccione (2002) for a slightly different form of coupling Hamiltonian but it is also applicable here in Anderson-Holstein model. It can also be checked that properties proposed in Breuer and Petruccione (2002) still hold:

Properties IV.1**.**

1. (I)

$\hat{D}(\omega)$ and $\hat{D}^{\dagger}(\omega)$ are eigen-operators of $\hat{H}_{s}$, namely,

[TABLE] 2. (II)

In the interaction picture, $\hat{D}_{I}(\omega,t)$ and $\hat{D}^{\dagger}_{I}(\omega,t)$ has the form

[TABLE] 3. (III)

$\hat{d}$ and $\hat{d}^{\dagger}$ can be decomposed into $\hat{D}(\omega)$ and $\hat{D}^{\dagger}(\omega)$ respectively. More specifically,

[TABLE] 4. (IV)

Then we can decompose coupling Hamiltonian $\hat{H}_{c}$ as

[TABLE]

and in the interaction picture

[TABLE]

These results directly follow from definition of $\hat{D}(\omega)$. Equation (IV.8) is essential in decomposing $\hat{d}(\hat{d}^{\dagger})$ in terms of levels of jumping. The reason that $\omega$ is the level of jumping can be observed from the definition that $\hat{D}(\omega)$ maps quantum state $\lvert\phi_{k^{\prime}}^{1},1\rangle$ to quantum state $\lvert\phi_{k}^{0},0\rangle$ where $k=k^{\prime}-\omega$. To prove the decomposition of $\hat{d}$ in terms of $\hat{D}(\omega)$, we use the completion relation $\sum_{k}(\Pi^{(0)}_{k}+\Pi^{(1)}_{k})=\text{Id}$,

[TABLE]

In the second step, we have used $\hat{d}\lvert\phi_{k^{\prime}}^{0},0\rangle=\hat{d}^{\dagger}\lvert\phi_{k}^{1},1\rangle=0$ for any $k,k^{\prime}\in\mathbb{N}$. It follows after taking the Hermitian conjugate that $\langle\phi_{k}^{1},1\rvert\hat{d}=0$; that is why there is only one term $\Pi^{(0)}_{k}\ \hat{d}\ \Pi^{(1)}_{k^{\prime}}$ left. In the third step, re-order the double summation is employed to first sum over all differences of levels, namely, $\omega$ and then sum over all possible combination of $k^{\prime},k\in\mathbb{N}$ such that $k^{\prime}-k=\omega$, where the latter sum gives $\hat{D}(\omega)$.

Starting from Equation (III.4), replacing $\hat{H}_{c,I}$ by Equation (IV.10) and opening the double commutators, we arrive at an alternative representation of Redfield equation

[TABLE]

where

[TABLE]

These two equations can be viewed as the Laplace transform of time correlation functions with frequency parameter $i\omega$.

IV.2 Secular approximation

From Equation (IV.11) and (IV.12), we observe that $\frac{d}{dt}\hat{\rho}_{s,I}(t)=\mathcal{O}\left(\frac{\alpha^{2}}{\epsilon}\right)$. This motivates the choice of the relaxation time as $\tau_{R}\equiv\frac{\epsilon}{\alpha^{2}}$. Recall that we have assumed $\alpha\ll\epsilon$ and $\alpha\ll 1$, then $\alpha^{2}\ll\epsilon$, or equivalently $\tau_{R}=\frac{\epsilon}{\alpha^{2}}\gg 1$. Take the integral of $\hat{\rho}_{s,I}(t)$ over time period $[t,t+r\tau_{R}]$ for $r=\mathcal{O}(1)$, we obtain

[TABLE]

We use the short-hand $Op(s)$ for simplicity to denote the long term involving operators in Equation (IV.11). Then change the variable $s=t+\tau_{R}s^{\prime}$,

[TABLE]

By Riemann-Lebesgue lemma, if $\omega^{\prime}\neq\omega$, $(\omega^{\prime}-\omega)\tau_{R}=\mathcal{O}(\tau_{R})\gg 1$,

[TABLE]

Then on the right hand side, terms involving $\omega^{\prime}-\omega\neq 0$ have negligible integral value. Hence,

[TABLE]

i.e.,

[TABLE]

This is known as secular approximation Breuer and Petruccione (2002), which we have justified here in the sense of coarse-grained approximation over relaxation time. Divide both side by $r\tau_{R}$ and then take the limit $r\rightarrow 0$, by fundamental theorem of calculus,

[TABLE]

Dropping the approximation, we arrived at the secular approximation, which is the basis for Lindblad equation:

[TABLE]

Remark*.*

By checking the previous argument, in fact, secular approximation is valid when $\frac{\alpha^{2}}{\epsilon}\ll 1$, that is, $\alpha^{2}\ll\epsilon$, which appears to be a weaker condition than $\alpha\ll\epsilon$ used for Born-Markov approximation.

IV.3 Lindblad equation in interaction picture

To write Equation (IV.14) in a Lindbladian form, we need to decompose coefficients $F(\omega)$ and $G(\omega)$ into their real and imaginary parts. Let $a_{F}(\omega):=F(\omega)+F(\omega)^{*}$ and $b_{F}(\omega):=\frac{F(\omega)-F(\omega)^{*}}{2i}$, then $F(\omega)\equiv\frac{a_{F}(\omega)}{2}+ib_{F}(\omega)$. Similarly, we can decompose $G(\omega)=\frac{\alpha_{G}(\omega)}{2}+ib_{G}(\omega)$. With these notations, Equation (IV.14) becomes Lindblad equation

[TABLE]

where Lindbladian correction Hamiltonian $\hat{\mathcal{H}}_{I}$ has the form

[TABLE]

and dissipative operator $\mathcal{D}$ has the form

[TABLE]

Recall that the general dissipative operator in Lindblad equation is a linear combination of

[TABLE]

where $\gamma$ is a constant Lindblad (1976). In Anderson-Holstein model, when $\gamma=a_{F}(\omega)$ the corresponding $\hat{L}=\hat{D}^{\dagger}(\omega)$; when $\gamma=a_{G}(\omega)$, the corresponding $\hat{L}=\hat{D}(\omega)$.

IV.4 Lindblad equation in Schrödinger picture

Transforming back into Schrödinger picture by $\hat{\rho}_{s}(t)=e^{-\frac{i}{\epsilon}\hat{H}_{s}t}\hat{\rho}_{s,I}(t)e^{\frac{i}{\epsilon}\hat{H}_{s}t}$, we obtain Lindblad equation in Schrödinger picture,

[TABLE]

where

[TABLE]

and

[TABLE]

Note that $\hat{\mathcal{H}}$ is invariant in different pictures since $\hat{D}(\omega)$ and $\hat{D}^{\dagger}(\omega)$ both appear in the same term and thus the factor $e^{\pm i\omega t}$ will always cancel during picture transformation; this cancellation also applies to the dissipative operator $\mathcal{D}$.

To understand the Lindbladian corrected Hamiltonian, we note that after some simple computation, it could be shown that

[TABLE]

Hence, for the new Hamiltonian, i.e., $\hat{H}_{s}+\alpha^{2}\hat{\mathcal{H}}$, the set of eigenstates are the same, $\lvert\phi_{k}^{0},0\rangle$ and $\lvert\phi_{k}^{1},1\rangle$ for $k\in\mathbb{N}$. In the new Hamiltonian, the energy eigenvalues, however, are perturbed by order $\mathcal{O}(\alpha^{2})$.

Physically, this perturbation of energy comes from the interaction between two quantum states $\lvert 0\rangle$ and $\lvert 1\rangle$ through the coupling with bath. More specifically, since the quantum eigenstate $\lvert\phi_{k}^{0},0\rangle$ interacts with quantum eigenstates $\lvert\phi_{k+\omega}^{1},1\rangle$ (for all possible $\omega\in\mathbb{Z}$), their interaction contributes to the change of energy; that is why there is term $\left|\langle\phi_{k}^{0}|\phi_{k+\omega}^{1}\rangle\right|^{2}$. The interaction is realized through the bath, hence the perturbed Hamiltonian should be weighted by time correlation functions, namely, terms $b_{F(G)}(\omega)$ and is also proportional to $\alpha^{2}$, the square of the coupling parameter.

The effect of Lindblad operator will be further investigated in Section VI below in the context of perturbation theory and it will be shown that Lindblad operator characterizes the hopping between quantum states $\lvert 0\rangle$ and $\lvert 1\rangle$. More specifically, the hopping rate out of eigenstate $\lvert\phi_{k}^{0},0\rangle$ is $\frac{\alpha^{2}}{\epsilon}\sum_{\omega}a_{F}(\omega)\left|\left\langle\phi_{k}^{0}|\phi_{k+\omega}^{1}\right\rangle\right|^{2}$. It is worth pointing out that this expression is quite similar to the perturbed energy eigenvalue as above, $\alpha^{2}\sum_{\omega}b_{F}(\omega)\left|\langle\phi_{k}^{0}|\phi_{k+\omega}^{1}\rangle\right|^{2}$. For the Laplace transform of time correlation functions, namely, $F(\omega)$ and $G(\omega)$ in Equation (IV.12), the real part contributes to (weak) hopping and the imaginary part contributes to (weak) perturbation to the energy eigenvalues.

IV.5 Coefficients and wide band approximation

It remains to determine the coefficients $a_{F,G}(\omega)$ and $b_{F,G}(\omega)$ in the Lindblad equation. Using oscillatory integral

[TABLE]

we could obtain that

[TABLE]

Thus matching the definition of $a_{F,G}(\omega)$ and $b_{F,G}(\omega)$’s

[TABLE]

Notice that $a_{F,G}(\omega)$ and $b_{F,G}(\omega)$ are (generalized) functions with respect to $\omega$. Even though we only need values at $\omega\in\mathbb{Z}$, but these functions are indeed well-defined on $\mathbb{R}$.

In Anderson-Holstein model, we have assumed that the electron bath is infinitely large, so continuum approximation appears to be a possible approach to simplify the coefficients. Assume that $V^{2}_{k}\equiv V^{2}(E_{k})$ is a continuous function of $E_{k}$. In the discrete case, suppose the total number of states in the bath is $N$, then $V^{2}(E)$ should be inversely proportional to $N$, to make the overall interaction strength between the system and bath remain at $\mathcal{O}(1)$: let $V^{2}(E)=\frac{\check{V}^{2}(E)}{N}$. Let $D$ be the energy band width of electron bath and $\nu(E)$ be the density of states at energy level $E$. In wide band approximation, to simplify the last equation (IV.21), it is assumed that the contribution to interaction strength from different energy levels of electron bath is approximately the same; explicitly, assume

[TABLE]

where $\Gamma$ is a constant Bruch et al. (2016). Then for a test function $g(\omega)$,

[TABLE]

Therefore, in the continuum limit,

[TABLE]

There are two ways to get rid of the characteristic function: the first way is to assume that $D=\infty$, mentioned in Ouyang et al. (2015); the second way is to consider $\epsilon\downarrow 0$. In either way, we end up with the approximation

[TABLE]

These two conditions are consistent with where the term $a_{F}(\omega)$ comes. $a_{F}(\omega)$ is part of interaction strength, which involves both electron bath and open quantum system; when the elctron bath is infinitely wide or the open quantum system falls into the semi-classical region, the (generalized) function $a_{F}(\omega)$ can be approximated in this way. Similarly, we can approximate

[TABLE]

For $b_{F}(\omega)$ and $b_{G}(\omega)$, we have not found easy expression for them.

V Semi-classical Limit

The semi-classical limit of Redfield equation has been proposed and studied in paper Dou et al. (2015a); Elste et al. (2008); the system of phase space functions obtained in the semiclassical limit of Redfield equation by applying Wigner transformation is called classical master equation. In the first part, we attempt to justify the formal derivation of Dou et al. (2015a) in a more mathematical way. In the second part, more importantly, we attempt to study the phase space counterparts of Lindblad equation by applying Wigner transformation. We call the system of phase space functions Lindbladian classical master equation. As far as we know, the Lindblad equation and its semi-classical limit for Anderson-Holstein model have not been studied.

V.1 Wigner transform, phase space functions and some notations

Recall that we have identified $\epsilon\equiv\frac{\hbar\omega_{s}}{\mathsf{E}}$ as semi-classical parameter. Also recall that the Wigner transformation of an operator $\hat{A}$ on $L^{2}(\mathbb{R})$ is defined by Wigner (1932); Zworski (2012)

[TABLE]

The subscript $\mathcal{W}$ indicates the Wigner transform. If $\hat{A}=\hat{\rho}$ is a density operator, it could be easily shown that $\int dxdp\ \frac{1}{2\pi\epsilon}\hat{\rho}_{\mathcal{W}}(x,p)=1$. Hence if we are applying Wigner transformation to density operators $\hat{\rho}$, the factor $\frac{1}{2\pi\epsilon}$ is needed to have integral equal to $1$ for function $\hat{\rho}_{\mathcal{W}}$.

For a single level quantum system, the phase space function for a quantum master equation is clear (see, for instance, Wigner (1932)). As for multi-level open quantum system, the definition for phase space functions is not very straightforward and we need to clarify this concept used below. For a general two-level open quantum system, the reduced density matrix can be written in the matrix form as

[TABLE]

For a general quantum master equation, it is expected that $\hat{\rho}_{0,1}(t)$ and $\hat{\rho}_{1,0}(t)$ do not vanish. However, it could be verified that if at time $t_{0}$, $\hat{\rho}_{s}(t)$ only have diagonal terms, i.e., $\hat{\rho}_{0,1}(t_{0})=\hat{\rho}_{1,0}(t_{0})=0$, then $\frac{d}{dt}\hat{\rho}_{0,1}(t_{0})=\frac{d}{dt}\hat{\rho}_{1,0}(t_{0})=0$ for both Redfield equation and Lindblad equation. Thus, $\hat{\rho}_{0,1}(t)=\hat{\rho}_{1,0}(t)=0$ for all $t$.

Therefore, in the below, we shall only consider diagonal elements of $\hat{\rho}_{s}(t)$, i.e., assume $\hat{\rho}_{s}(t)=\hat{\rho}_{0}(t)\lvert 0\rangle\langle 0\rvert+\hat{\rho}_{1}(t)\lvert 1\rangle\langle 1\rvert$ for all $t$. In matrix form,

[TABLE]

The phase space function by applying Wigner transformation for $\hat{\rho}_{m}(t)$ is denoted by $\varrho_{m}(x,p,t)$, for $m\in\{0,1\}$, namely,

[TABLE]

Hence, it is not difficult to show that $\int dxdp\ \varrho_{0}(x,p,t)+\varrho_{1}(x,p,t)=1$.

We need to use the following lemma below.

Lemma 1**.**

(Semi-classical expansion, (Zworski, 2012, p. 66 - 68))

(\hat{A}\hat{B})_{\mathcal{W}}=\hat{A}_{\mathcal{W}}\hat{B}_{\mathcal{W}}+\frac{i\epsilon}{2}\mathopen{\big{missing}}\{\hat{A}_{\mathcal{W}},\hat{B}_{\mathcal{W}}\mathopen{\big{missing}}\}+\mathcal{O}(\epsilon^{2})* where $\hat{A}$ and $\hat{B}$ are two operators. $(\hat{A}\hat{B})_{\mathcal{W}}$ is linear with respect to both $\hat{A}_{\mathcal{W}}$ and $\hat{B}_{\mathcal{W}}$.* 2. 2.

In particular, \mathopen{\big{missing}}[\hat{A},\hat{B}\mathopen{\big{missing}}]_{\mathcal{W}}=i\epsilon\mathopen{\big{missing}}\{\hat{A}_{\mathcal{W}},\hat{B}_{\mathcal{W}}\mathopen{\big{missing}}\}+\mathcal{O}(\epsilon^{3}) where \mathopen{\big{missing}}\{\cdot,\cdot\mathopen{\big{missing}}\} is the Poisson bracket, defined by111We choose to follow the usual convention here, which differs by a negative sign compared to Zworski (2012)**

[TABLE]

V.2 Classical master equation

As discussed in last Subsection, we shall only consider the reduced density matrix of the form $\hat{\rho}_{s}(t)=\hat{\rho}_{0}(t)\lvert 0\rangle\langle 0\rvert+\hat{\rho}_{1}(t)\lvert 1\rangle\langle 1\rvert$. From Equation (III.7), we could derive that the time evolution equations for $\hat{\rho}_{0}(t)$ and $\hat{\rho}_{1}(t)$ is

[TABLE]

and

[TABLE]

which agree with Equation (14) and (15) in Ref Dou et al. (2015c).

Using Lemma 1, we can calculate the equations for the corresponding Wigner transformation as

[TABLE]

For clarity in equation, the coordinates $(x,p)$ are omitted in phase space functions and in $H_{0}$, $H_{1}$ as well. After dropping higher order terms of $\mathcal{O}(\alpha^{2})$

[TABLE]

where hopping rates

[TABLE]

Notice that hopping rates are functions of phase space coordinates $x$ and $p$. They describe how fast the jumping between states $\lvert 0\rangle$ and $\lvert 1\rangle$ depending on $(x,p)$.

Remark*.*

Compared with the result in Dou et al. (2015a), the rates $\gamma_{0\to 1}$ and $\gamma_{1\to 0}$ we have above are considerably more complicated. Here we provide a heuristic simplification of the expression, though we do not know how to justify the argument on a more rigorous level.

Equation (60) in Curtright et al. (2014) shows that if $\hat{H}$ is a Hamiltonian for harmonic oscillators,

[TABLE]

where $H=(\hat{H})_{\mathcal{W}}$. Since $\hat{H}_{0}$ and $\hat{H}_{1}$ are Hamiltonian for harmonic oscillators, by using the last equation,

[TABLE]

In the last step, we use $\cos(\tau/2)\rightarrow 1$ and $\tan(\tau/2)\rightarrow\tau/2$ as $\tau\rightarrow 0$. Then, hopping rates can be written as

[TABLE]

Heuristically, if we approximate $(e^{-i\hat{H}_{1}\tau})_{\mathcal{W}}(e^{i\hat{H}_{0}\tau})_{\mathcal{W}}$ by $e^{-iU(x)\tau}$, i.e., assume that the integral is mostly contributed from $\tau$ near [math], then

[TABLE]

which becomes the result in Dou et al. (2015a). $U^{-1}(x)$ is the inverse function of $U(x)$; since $U(x)$ is a linear function, $U^{-1}(x)$ is well-defined. The second and third step of last equation use similar computation as wide band approximation, which has been shown in details in Section IV.5. This heuristic computation gives a nice simple expression for the rates, but it should be pointed out we do not know how to justify the crucial approximation above of $(e^{-i\hat{H}_{1}\tau})_{\mathcal{W}}(e^{i\hat{H}_{0}\tau})_{\mathcal{W}}$ by $e^{-iU(x)\tau}$.

V.3 Lindbladian classical master equation

Recall that if we assume at time $t_{0}$, $\hat{\rho}_{s}(t_{0})$ is diagonal, then reduced density operator $\hat{\rho}_{s}(t)=\hat{\rho}_{0}(t)\lvert 0\rangle\langle 0\rvert+\hat{\rho}_{1}(t)\lvert 1\rangle\langle 1\rvert$ without terms involving $\lvert 1\rangle\langle 0\rvert$ nor $\lvert 0\rangle\langle 1\rvert$. The time-evolution Lindblad equations can be written more explicitly as

[TABLE]

[TABLE]

The system of time-evolution equations obtained by applying Wigner transformation to Lindblad equation is given by, after some straightforward calculations

[TABLE]

[TABLE]

where

[TABLE]

The error terms with order $\mathcal{O}(\alpha^{2})$ come from the Wigner transformation of hopping terms and the error terms with order $\mathcal{O}(\alpha^{2}\epsilon^{2})$ come from Wigner transformation of commutators.

VI Comparison of Redfield equation and Lindblad equation from perturbation theory

In this section, we intend to use perturbation theory to understand the similarity and difference of Redfield equation and Lindblad equation, by considering the hopping between different quantum states. The mathematical tool for perturbation theory has been widely studied, and one of the most famous formulas from that is Fermi Golden Rule.

The main finding is that if the quantum system is prepared at a pure state $\lvert\phi_{i}^{0},0\rangle$, then its hopping rates to other quantum states $\lvert\phi_{f}^{0},0\rangle$ and $\lvert\phi_{f}^{1},1\rangle$ are the same up to the first order, for Redfield and Lindblad equations. By linearity of Redfield and Lindblad equation, this conclusion is also true if initially $\hat{\rho}_{m}$ ($m=0,1$) are diagonal.

VI.1 Perturbation of Redfield and Lindblad equations

As discussed above, Lindblad equation can be derived by using secular approximation from Redfield equation, under the condition that $\alpha^{2}\ll\epsilon$. Under this condition, the hopping term in Redfield equation (Equation (IV.11)) could be considered as a small perturbation, that is, $\frac{\alpha^{2}}{\epsilon}$ is regarded as perturbation parameter below.

In Schrödinger picture, the Redfield equation has the form, by using Equation (IV.11),

[TABLE]

After rewriting it in terms of $a_{F(G)}(\omega)$ and $b_{F(G)}(\omega)$, it becomes

[TABLE]

where operator $\mathcal{R}$ has the form

[TABLE]

If we only keep those terms with $\omega^{\prime}=\omega$, it becomes Lindblad equation, that is Equation (IV.18).

Then, consider the perturbation for the system initially at $\lvert\phi_{i}^{0},0\rangle$ and ending at $\lvert\phi_{f}^{0},0\rangle$ or $\lvert\phi_{f}^{1},1\rangle$. We are interested in finding the transition between different quantum states. Because of the symmetry of this problem, it suffices to look at transition from $\lvert\phi_{i}^{0},0\rangle$ to $\lvert\phi_{f}^{0},0\rangle$ and to $\lvert\phi_{f}^{1},1\rangle$.

Assume that

[TABLE]

Then by collecting terms of the same order,

[TABLE]

For the leading order, the solution is trivial to find $\hat{\rho}_{s,0}(t)=\lvert\phi_{i}^{0},0\rangle\langle\phi_{i}^{0},0\rvert\equiv\Pi_{i}^{(0)}$. For the first order, after some computing,

[TABLE]

The terms involving $a_{G}(\omega)$ and $b_{G}(\omega)$ vanish, since there is no hopping from state $\lvert 1\rangle$ to $\lvert 0\rangle$ in our setting up. To study transition between different quantum states, let

[TABLE]

where $m=0,1,2,\cdots$ represents the perturbation order. For the leading order, apparently, $\lambda_{f}^{(0)}(t)=\delta_{i,f}$ and $\theta^{(0)}_{f}(t)=0$ for $f\in\mathbb{N}$.

As for the first-order, it could be obtained that

[TABLE]

and

[TABLE]

As we observe, for the diagonal elements of $\mathcal{R}(\hat{\rho}_{s,0}(t))$, the only contribution comes from terms with $\omega=\omega^{\prime}$. Therefore, if this condition $\omega^{\prime}=\omega$ is imposed beforehand, the result for $\frac{d}{dt}\lambda_{f}^{(1)}(t)$ and $\frac{d}{dt}\theta_{f}^{(1)}(t)$ will not change. This implies that if we are applying perturbation to Lindblad equation, we should obtain exactly the same set of equations for the diagonal elements up to the first order.

It could be noticed that the rate of probability leaving $\lvert\phi_{i}^{0},0\rangle$ is the summation of rate of probability entering $\lvert\phi_{f}^{1},1\rangle$ for varying $f\in\mathbb{N}$. The rate of leaving $\lvert\phi_{i}^{0},0\rangle$ is $\frac{\alpha^{2}}{\epsilon}\sum_{\omega}a_{F}(\omega)\left|\langle\phi_{i}^{0}\lvert\phi_{i+\omega}^{1}\rangle\right|^{2}$ and it hops to state $\lvert\phi_{f}^{1},1\rangle$ with rate $\frac{\alpha^{2}}{\epsilon}a_{F}(\omega)\left|\langle\phi_{i}^{0}\lvert\phi_{i+\omega}^{1}\rangle\right|^{2}$ if $f=i+\omega$ for some $\omega\in\mathbb{Z}$.

VI.2 More general initial condition

A slightly more general initial condition is to assume that $\hat{\rho}_{0}$ and $\hat{\rho}_{1}$ are diagonal, i.e., $\hat{\rho}_{0}(0)=\sum_{k}\lambda_{k}\lvert\phi_{k}^{0}\rangle\langle\phi_{k}^{0}\rvert$ and $\hat{\rho}_{1}(0)=\sum_{k}\theta_{k}\lvert\phi_{k}^{1}\rangle\langle\phi_{k}^{1}\rvert$. Then the leading order is

[TABLE]

Then we need to compute the first order,

[TABLE]

Then, after some computation

[TABLE]

and

[TABLE]

Hence the hopping rate from $\lvert\phi_{i}^{0},0\rangle$ to $\lvert\phi_{f}^{1},1\rangle$ is $\frac{\alpha^{2}}{\epsilon}a_{F}(f-i)\left|\langle\phi_{i}^{0}|\phi_{f}^{1}\rangle\right|^{2}$ and hopping rate from $\lvert\phi_{i}^{1},1\rangle$ to $\lvert\phi_{f}^{0},0\rangle$ is $\frac{\alpha^{2}}{\epsilon}a_{G}\left(i-f\right)\left|\langle\phi_{f}^{0}|\phi_{i}^{1}\rangle\right|^{2}$. Because $\lambda_{k}(t)$ and $\theta_{k}(t)$ have an interpretation as the probability at state $\lvert\phi_{k}^{0},0\rangle$ or $\lvert\phi_{k}^{1},1\rangle$, then Equations (VI.7) and (VI.8) can be interpreted as the Kolmogorov’s backward equation for a continuous time Markov Chain. To make sure that the first-order perturbation is valid, the time horizon is assumed to be $t\ll\frac{\epsilon}{\alpha^{2}}$.

The computation shows that if initially $\hat{\rho}_{0}(0)$ and $\hat{\rho}_{1}(0)$ are diagonal, then Redfield and Lindblad equation lead to the same transition rate up to the first order. Since the operator $\mathcal{R}$ is linear with respect to the input density operator, this result is just a natural extension of that in last Subsection.

VI.3 Interpretation of hopping rate in LCME

We continue to assume that $\hat{\rho}_{m}(0)$ ($m=0,1$) are diagonal and in this Subsection, we will consider Lindblad equation only. It could be easily verified that $\hat{\rho}_{m}(t)$ are diagonal for any time $t\geq 0$. Hence

[TABLE]

Consider the jumping leaving $\lvert 0\rangle$ to $\lvert 1\rangle$ in LCME, that is,

[TABLE]

in Equation (V.7). By using definition of $\hat{D}^{(\dagger)}(\omega)$ in Equation (IV.4), and using Lemma 1, we could obtain

[TABLE]

This equation shows that the hopping rate from state $\lvert 0\rangle$ to $\lvert 1\rangle$ in semi-classical limit is the summation of the contribution from each state $\lvert\phi_{k}^{0},0\rangle$. More specifically, the probability at state $\lvert\phi_{k}^{0},0\rangle$ in the semi-classical sense is $\lambda_{k}(t)\frac{1}{2\pi\epsilon}\left(\lvert\phi_{k}^{0}\rangle\langle\phi_{k}^{0}\rvert\right)_{\mathcal{W}}$; the hopping rate out of state $\lvert\phi_{k}^{0}\rangle$ is $\left(\frac{\alpha^{2}}{\epsilon}\sum_{\omega}a_{F}(\omega)\left|\langle\phi_{k}^{0}\rvert\phi^{1}_{k+\omega}\rangle\right|^{2}\right)$; their product is exactly the contribution from quantum state $\lvert\phi_{k}^{0},0\rangle$.

Notice that the hopping rate obtained here is consistent with Equation (VI.5) obtained by perturbation theory. This matches our intuition and it connects the hopping rate in LCME with the hopping rate from Lindblad equation.

VI.4 Discussion on Franck-Condon Blockade

Using the wide band approximation, the hopping rate could be explicitly computed using Franck-Condon factors Koch and von Oppen (2005); Dou et al. (2015c). The Franck-Condon factor is

[TABLE]

where $L_{k}^{(k^{\prime})}(x)$ is the generalized Laguerre polynomial, $N=\min(n,m)$ and $M=\max(n,m)$. Then the hopping rate from $\lvert\phi_{n}^{0},0\rangle$ to $\lvert\phi_{m}^{1},1\rangle$ is

[TABLE]

We observe that the hopping rate is small due to the weak phonon-reservoir coupling (i.e., $\frac{\alpha^{2}}{\epsilon}\ll 1$) and also since the Franck-Condon factor $\exp(-g^{2}/\epsilon)$ becomes exponentially small as $g^{2}/\epsilon\rightarrow\infty$, known as the Franck-Condon blockade Koch and von Oppen (2005); Ryndyk (2016) for large on-the-molecule electron-phonon coupling (i.e., $g\gg 1$). From the last expression, we notice that the Franck-Condon blockade occurs when the ratio of electron-phonon coupling rate $g^{2}$ and semi-classical parameter $\epsilon$ is large; in particular, even for finite $g$, if $\epsilon\rightarrow 0$, such exponentially small hopping rate also appears.

VII Discussion and Conclusion

In this paper, we have revisited the derivation of Redfield equation by using time-convolutionless equation (TCL) and consequently derive the Lindblad equation using secular approximation, in the context of Anderson-Holstein model. As an analogy to classical master equation, Lindbladian classical master equation (LCME) is introduced and its form is given in Equations (V.7) and (V.8). The comparison between Redfield equation and Lindblad equation from the perspective of perturbation theory is considered: Redfield equation and Lindblad equation yield the same hopping rate in the first order; in other words, the dynamics of their diagonal elements of reduced density operator $\hat{\rho}_{s}$ are the same, to the first order. The condition of the derivation and perturbation result both suggest that Lindblad equation might be a better candidate for studying Anderson-Holstein model than Redfield equation and reasons are listed as follow: first, from the derivation, they are at the same level, i.e., both under weak-coupling limit $\alpha\ll\epsilon$; second, the condition of deriving Lindblad equation from Redfield equation is redundant and there is, in fact, no further constraint in approximating Redfield equation by Lindblad equation; third, they have the same hopping rates between eigenstates up to the first order in perturbation theory; finally, the analysis for Lindblad equation is easier.

Usually, for multi-level open quantum systems, the density operator $\hat{\rho}_{s}(t)$ does not necessarily have vanishing $\hat{\rho}_{0,1}(t)$ and $\hat{\rho}_{1,0}(t)$. These two off-diagonal operators will cause more challenge in both analysis and numerics. As we have shown, for Redfield equation, if we start from diagonal density operator, i.e., $\hat{\rho}_{s}(t)=\hat{\rho}_{0}(t)\lvert 0\rangle\langle 0\rvert+\hat{\rho}_{1}(t)\lvert 1\rangle\langle 1\rvert$ at some time $t=t_{0}$, then it remains to be in this form for all time $t$. If we start from a general density operator $\hat{\rho}_{s}(t)$, then analyzing the time-evolution of off-diagonal terms (i.e., $\hat{\rho}_{0,1}(t)$ and $\hat{\rho}_{1,0}(t)$) would be challenging and it could be our next stage of research. For Lindblad equation, the same phenomenon appears. This approach of only considering diagonal operators with vanishing $\hat{\rho}_{0,1}(t)$ and $\hat{\rho}_{1,0}(t)$, has its restriction but it could lead into simple equations after applying Wigner transformation. The semi-classical limits for Redfield equation and Lindblad equation have similar form.

Some continuing works could include (1) how to perform further approximations to simplify the hopping coefficients in LCME to obtain a simpler equation; (2) a systematic and controllable numerical method to solve Lindblad equation in infinite dimensional Hilbert space L^{2}(\mathbb{R})\otimes\text{span}\mathopen{\big{missing}}\{\lvert 0\rangle,\lvert 1\rangle\mathclose{\big{missing}}\}. We will leave these for future works.