Electron Counting Statistics for Non-Additive Environments

TL;DR

This paper develops a model for electron transport in molecular electronics considering strong vibrational coupling, revealing non-additive effects on current, noise, and Fano factor using counting statistics.

Contribution

It introduces a non-additive approach to model lead couplings in systems with strong vibrational interactions, advancing understanding of electron transport in molecular devices.

Findings

Strong vibrational coupling affects electron flow statistics.

Non-additive lead couplings alter current and noise characteristics.

The model captures non-Markovian effects in electron transport.

Abstract

Molecular electronics is a rapidly developing field focused on using molecules as the structural basis for electronic components. It is common in such devices for the system of interest to couple simultaneously to multiple environments. Here we consider a model comprised of a double quantum dot (or molecule) coupled strongly to vibrations and weakly to two electronic leads held at arbitrary bias voltage. The strong vibrational coupling invalidates treating the bosonic and electronic environments simply as acting additively, as would be the case in the weak coupling regime or for flat leads at infinite bias. Instead, making use of the reaction coordinate framework we incorporate the dominant vibrational coupling effects within an enlarged system Hamiltonian. This allows us to derive a non-additive form for the lead couplings that accounts properly for the influence of strong and non-Markovian coupling between the double dot system and the vibrations. Applying counting statistics techniques we track electron flow between the double dot and the electronic leads, revealing both strong-coupling and non-additive effects in the electron current, noise and Fano factor.

Equations

$$H=H_{\mathrm{S}}+H_{{\mathrm{E}}_{\mathrm{ph}}}+H_{{\mathrm{E}}_{\mathrm{L}}}+H_{{\mathrm{E}}_{\mathrm{R}}}+H_{{\mathrm{I}}_{\mathrm{ph}}}+H_{{\mathrm{I}}_{\mathrm{L}}}+H_{{\mathrm{I}}_{\mathrm{R}}},$$

$$H_{\mathrm{S}}={ϵ}_{\mathrm{L}}|L⟩⟨L|+{ϵ}_{\mathrm{R}}|R⟩⟨R|+T\left(|L⟩⟨R|+|R⟩⟨L|\right)$$

$$H_{{\mathrm{E}}_{\mathrm{ph}}}$$

$$H_{{\mathrm{E}}_j}$$

$$H_{{\mathrm{I}}_{\mathrm{ph}}}$$

$$H_{{\mathrm{I}}_j}$$

$$|{{\rm Y}}_1⟩$$

$$|{{\rm Y}}_2⟩$$

$${{\rm Y}}_1$$

$${{\rm Y}}_2$$

$$H_{{\mathrm{I}}_{\mathrm{L}}}=A_1\otimes B_1+A_2\otimes B_2,$$

$${\mathcal{L}}_{\mathrm{L}}{\rho }_{\mathrm{S}}(t)$$

$$+[{\rho }_{\mathrm{S}}(t)A_1(-\tau ),A_2]C_{12}(-\tau )$$

$$+[A_2,A_1(-\tau ){\rho }_{\mathrm{S}}(t)]C_{21}(\tau )$$

$$+[{\rho }_{\mathrm{S}}(t)A_2(-\tau ),A_1]C_{21}(-\tau ),$$

$$A_1(-\tau )$$

$$A_2(-\tau )$$

$$C_{12}(\tau )$$

$$C_{21}(\tau )$$

$${\mathcal{L}}_{\mathrm{R}}{\rho }_{\mathrm{S}}(t)$$

$$+[{\rho }_{\mathrm{S}}(t)A_3(-\tau ),A_4]C_{34}(-\tau )$$

$$+[A_4,A_3(-\tau ){\rho }_{\mathrm{S}}(t)]C_{43}(\tau )$$

$$+[{\rho }_{\mathrm{S}}(t)A_4(-\tau ),A_3]C_{43}(-\tau ),$$

$$A_3(-\tau )$$

$$A_4(-\tau )$$

$$C_{34}(\tau )$$

$$C_{43}(\tau )$$

$$H_{{\mathrm{I}}_{\mathrm{ph}}}=A_5\otimes B_5,$$

$$A_5(t)$$

$$\displaystyle+\zeta T{\rm sin}(2\zeta t)\sigma_{y}\Big{]},$$

$$\displaystyle\mathcal{L}_{\mathrm{ph}}\rho_{\rm S}(t)=-\Big{\{}$$

$$\displaystyle+M_{x}[\sigma_{z},[\sigma_{x},\rho_{\rm S}(t)]]\Big{\}},$$

$$J(\omega )=\sum _q|g_q{|}^2\delta (\omega -{\omega }_q)=\frac{1}{\pi }\frac{2\lambda \omega {\omega }_0^2\gamma }{({\omega }_0^2-{\omega }^2)+{\gamma }^2{\omega }^2}.$$

$${\mathcal{J}}_{L/R}(\omega )=\sum _{k_{L/R}}|t_{k_{L/R}}{|}^2\delta (\omega -{ϵ}_{\mathrm{L}/R}),$$

$${\dot{\rho }}_{\mathrm{S}}(t)$$

$$H$$

$$\to H_{{\mathrm{S}}^{\prime }}+H_{{\mathrm{E}}^{\prime }}+H_{{\mathrm{I}}^{\prime }}+H_{{\mathrm{E}}_{\mathrm{L}}}+H_{{\mathrm{E}}_{\mathrm{R}}}+H_{{\mathrm{I}}_{\mathrm{L}}}+H_{{\mathrm{I}}_{\mathrm{R}}}.$$

$$H_{{\mathrm{S}}^{\prime }}$$

$$+\kappa (|L⟩⟨L|-|R⟩⟨R|)(a^{†}+a)+\Omega a^{†}a,$$

$$\kappa (a^{†}+a)=\sum _q{g}_q(a_q^{†}+a_q),$$

$$H_{{\mathrm{E}}^{\prime }}=\sum _k{\nu }_k{b}_k^{†}{b}_k,$$

Full text

I Introduction

The original idea for a molecular nanojunction was put forward by Aviram and Ratner in 1974 Aviram and Ratner (1974), in which they suggested using a single organic molecule to develop a simple electronic rectifier. Since this initial proposal the field of molecular electronics has grown quickly, both in the construction of electronic components from molecular building blocks, and in the modelling and measurement of devices to tailor and probe their useful characteristics Tao (2006); Sun et al. (2014); Aradhya and Venkataraman (2013); Chen et al. (2007); Nichols and Higgins (2015). The past few years have seen numerous breakthroughs in the production of molecular nanojunctions [Reed et al., 1997; Bezryadin et al., 1997; Park et al., 1999; Kubatkin et al., 2003; Cui et al., 2001; Nicewarner-Peña et al., 2001; Nichols and Higgins, 2015; Sun et al., 2014], while theoretical advances have also been made with recent works examining a wide variety of nanojunction attributes; examples include current-voltage Tao (2006), conductance Nitzan and Ratner (2003); Chen et al. (2007) and noise (Fano factor) Koch and von Oppen (2005a) characteristics, along with thermoelectric, optical and electromechanical processes Aradhya and Venkataraman (2013); Cui et al. (2017); Segal and Agarwalla (2016); Moghaddasi Fereidani and Segal (2019). One aspect of key importance in terms of further experimental progress and the viability of applications is understanding the role of electron-vibrational coupling Kushmerick et al. (2004); Koch and von Oppen (2005a, b); Galperin et al. (2006); Kießlich et al. (2007); Hübener and Brandes (2009); Härtle et al. (2009); Härtle and Thoss (2011); Emary et al. (2012); Santamore et al. (2013); Pigeon et al. (2017); Stones and Olaya-Castro (2017); Sowa et al. (2017a, b). This has been explored, for example, through nanojunction models with coupling to individual vibrational modes Koch and von Oppen (2005a); Harbola et al. (2006a); Santamore et al. (2013); Schaller (2014); You et al. (2017). Our focus here is on a continuum of modes, where we wish to avoid either a weak electron-vibrational coupling approximation or the infinite lead bias regime, both of which can be used to simplify the resulting formalism. We thus aim to present an analysis of nanojunction models, and specifically electron counting statistics techniques, within arbitrary bias regimes for strong and non-Markovian vibrational coupling Strasberg et al. (2016); Sowa et al. (2017b, a); Krause et al. (2015).

We shall consider a molecular double quantum dot (DQD) system comprised of two electronic levels, which are coupled to a pair of fermionic leads and a bosonic vibrational (phonon) environment. This set up is a working model for a molecular nanojunction in which vibrational influences are significant. We address the challenge of strong electron-phonon coupling using the reaction coordinate (RC) mapping Iles-Smith et al. (2014, 2016); Strasberg et al. (2016); Maguire et al. (2018). This is a (unitary) Hamiltonian level mapping that transforms the coupling of the DQD electronic system to the original continuum phonon environment to coupling only to a single collective mode (the RC), which in turn interacts weakly with a residual bosonic environment. It thus allows for the inclusion of the RC within an enlarged (augmented) system along with the DQD to capture the dominant strong-coupling and non-Markovian effects. We then treat the remaining (residual) bosonic environment and the fermionic leads as being weakly coupled to the augmented system using standard Born-Markov (BM) master equation methods.

As a means to compare with common approximations, we shall analyse our model within three approaches. We begin with a standard weak coupling theory, in which we treat the interactions between our system and each of the environments using the BM approximations without any RC mapping. We shall compare this weak coupling approach with two strong-coupling methods applied to the same underlying Hamiltonian. In the first of these we consider the phonons within the RC mapped formalism to incorporate strong electron-phonon interactions. However, the leads will be treated through standard dissipators having no indication of this strong coupling, i.e. they will be incorporated additively. In the second strong coupling approach we again consider the phonon bath within the RC formalism, however, unlike before, we shall then consistently derive the lead dissipators as being weakly coupled to our enlarged augmented system that includes both the DQD and the RC. As such the leads become “aware” of the phonon coupling within the augmented system, resulting in a non-additive treatment. In this way we can build up our understanding of the system from the weak coupling theory, to a strong coupling additive theory, and then finally arrive at a strong coupling theory containing non-additive contributions. This allows us to highlight and isolate both strong-coupling and non-additive vibrational contributions to the system transport properties.

We shall employ full counting statistics techniques to probe electron transport, which enables calculation of the zero frequency current, noise and Fano factor of our DQD system Flindt et al. (2005, 2010). From these quantities we can ascertain information regarding the electron correlations present within the probability distribution of our system. Originally, these techniques were used alongside perturbative treatments for the phonon environment Kießlich et al. (2007); Schaller (2014). However recent works have successfully applied electron counting statistics to systems where the electron-phonon interaction is non-perturbative, for example through the hierarchical equations of motion approach Stones and Olaya-Castro (2017); Cerrillo et al. (2016) or a polaron transformation Sowa et al. (2017a). One feature of the RC treatment worth noting is that although our electronic system dynamics is non-Markovian, we work with a Markovian master equation in the augmented (RC mapped) representation, and as such the electron counting statistics formalism is easily applied.

Molecular nanojunction models are also commonly probed within the infinite bias regime, whereby the leads are held at chemical potentials that ensure a unidirectional flow of electrons. For flat lead spectral densities the infinite bias case is necessarily additive (see also below), however the RC formalism allows us to relax this boundary, thus providing insight into finite bias regimes and the impact of non-additive lead-phonon contributions on the transport statistics.

The paper is organised as follows. In Section II we shall introduce the theoretical model, briefly commenting on each component of the Hamiltonian, before outlining the derivation of the weak-coupling treatment in Section II.1. In Section II.2 we shall discuss the implementation of the RC formalism for strong phonon coupling and derive the lead components for the non-additive master equation. Section II.3 describes the procedure for implementing full counting statistics within our system. In Section III we shall present our results. We demonstrate that non-additive contributions have an impact at finite bias when considering the zero frequency current, noise, and Fano factor for flat leads, and even also at infinite bias if the leads are considered to have Lorentzian profiles. Finally, in Section IV we shall summarise our findings and discuss potential further work.

II Theoretical Model

A schematic representation of our model system is given in Fig. 1. We consider a DQD that is coupled to a phonon environment and two fermionic leads. Each dot is represented by an electronic level with energies $\epsilon_{\rm L}$ and $\epsilon_{\rm R}$ for the left and right dot, respectively. Throughout this work we impose the Coulomb blockade limit, ensuring that only one electron is present on either dot at any given time. We thus work in a reduced system space spanned by $\{\ket{0},\ket{L},\ket{R}\}$, which we term the site basis, where $\ket{0}$ represents the case in which neither dot is occupied and $\ket{L}$ ($\ket{R}$) denotes an electron present on the left (right) dot. The dots are coupled by a tunneling parameter $T$, which provides the primary transport channel within our system, while the phonon bath provides a secondary channel of transport provided that the tunneling amplitude is non-zero.

After applying a Jordan-Wigner transformation Schaller (2014) and restricting to the Coulomb blockade regime, the Hamiltonian governing the full system is given by (we take $\hbar=1$ throughout)

[TABLE]

where

[TABLE]

represents the DQD system. The Hamiltonian governing the environments has been split into components $H_{\rm E_{ph}}$ representing the bosonic phonon bath contributions and $H_{{\rm E}_{j}}$ ($j={\rm L,R}$) representing the fermionic leads:

[TABLE]

We consider phonon modes $q$ of frequencies $\omega_{q}$ and standard bosonic creation (annihilation) operators $a_{q}^{\dagger}$ ($a_{q}$). Similarly, the lead modes $k_{j}$ correspond to energy levels $\epsilon_{k_{j}}$ and standard fermionic creation (annihilation) operators $c_{k_{j}}^{\dagger}$ ($c_{k_{j}}$). Finally, the interaction Hamiltonians are given explicitly as (again $j={\rm L,R}$)

[TABLE]

Here, each phonon mode $q$ is assumed to couple with strength $g_{q}$ to the dot occupation difference Kießlich et al. (2007); Stones and Olaya-Castro (2017). The tunnel coupling strengths between each dot $j$ and the associated lead modes $k_{j}$ are labelled $t_{k_{j}}$.

II.1 Weak-coupling master equation

We begin by presenting a standard BM master equation for the DQD, whereby all environments are considered to be weakly interacting with the system. This will form one aspect of our comparisons later on, in order to highlight the impact of the strong phonon coupling physics captured within the RC treatment. We derive each environmental contribution individually in this case, as they are additive within the BM approximations. Beginning with the lead contributions, the nanojunction has two tunneling transport channels defined through the energy eigenstates of our system Hamiltonian, $\ket{\Upsilon_{1}}$ and $\ket{\Upsilon_{2}}$, which have energies $\Upsilon_{1}$ and $\Upsilon_{2}$, respectively:

[TABLE]

Here, $\Delta=({\epsilon_{\rm L}-\epsilon_{\rm R}})/{2}$ and $\alpha=\Delta-\sqrt{\Delta^{2}+T^{2}}$. In order to derive a standard BM master equation we express our site basis states $\ket{L}$ and $\ket{R}$ in terms of these eigenstates. The left lead interaction Hamiltonian can be written as

[TABLE]

where $A_{1}=\alpha\ket{0}\bra{\Upsilon_{2}}-T\ket{0}\bra{\Upsilon_{1}}$, $A_{2}=\alpha\ket{\Upsilon_{2}}\bra{0}-T\ket{\Upsilon_{1}}\bra{0}$, $B_{1}=\sum_{k_{\rm L}}t_{k_{\rm L}}c^{\dagger}_{k_{\rm L}}$ and $B_{2}=\sum_{k_{\rm L}}t^{*}_{k_{\rm L}}c_{k_{\rm L}}$. Moving into the interaction picture with respect to $H_{\rm S}+H_{\rm E_{\rm L}}$ we trace out the lead degrees of freedom within the standard BM approximations. On returning to the Schrödinger picture we find the left lead weak-coupling Liouvillian

[TABLE]

where $\rho_{\rm S}(t)$ is the reduced density operator of the DQD system. Here, interaction picture system operators and lead correlation functions are given as

[TABLE]

with $f_{j}(\epsilon)=(e^{\beta_{\rm el}(\epsilon-\mu_{j})}+1)^{-1}$ for lead $j$ with chemical potential $\mu_{j}$ and inverse temperature $\beta_{\rm el}$. Further details of this procedure are given in Appendix A, as is the explicit form of $\mathcal{L}_{\mathrm{L}}\rho_{\rm S}(t)$ in the site basis.

For the right lead contributions we have an analogous result

[TABLE]

with

[TABLE]

Again, the explicit form of $\mathcal{L}_{\mathrm{R}}\rho_{\rm S}(t)$ in the site basis is given in Appendix A.

The phonon contributions are more easily evaluated directly in terms of the site basis states, where

[TABLE]

with $A_{5}=\ket{L}\bra{L}-\ket{R}\bra{R}$ and $B_{5}=\sum_{q}g_{q}(a^{\dagger}_{q}+a_{q})$. Moving into the interaction picture with respect to $H_{\rm S}+H_{\rm E_{ph}}$ and making use of the Baker-Campbell-Hausdorff formula Sattinger and Weaver (2013), we find

[TABLE]

where $\zeta=\sqrt{\Delta^{2}+T^{2}}$ and $\sigma_{i}$ are the usual Pauli matrices in the system subspace $\{\ket{L},\ket{R}\}$, with $\sigma_{z}=\ket{L}\bra{L}-\ket{R}\bra{R}$. Following the standard BM procedure once more, we derive the weak-coupling phonon Liouvillian

[TABLE]

where $M_{z}=\frac{2\Delta^{2}\pi\lambda}{\zeta^{2}\beta_{\rm ph}\omega_{0}}+\frac{2\pi T^{2}}{\zeta^{2}}J(2\zeta)\mathrm{coth}(\beta_{\rm ph}\zeta)$, $M_{y}=\frac{2i\pi TJ(2\zeta)}{\zeta}$ and $M_{x}=\frac{2\Delta T\pi\lambda}{\zeta^{2}\beta_{\rm ph}\omega_{0}}-\frac{2\pi\Delta T}{\zeta^{2}}J(2\zeta)\mathrm{coth}(\beta_{\rm ph}E)$, for inverse phonon temperature $\beta_{\rm ph}$. Throughout this work we consider a Drude-Lorentz spectral density for the phonon bath, defined as

[TABLE]

Here $\lambda$ is the reorganisation energy of the phonons, $\omega_{0}$ is the central peak and $\gamma$ is a parameter which tunes the spectral density from sharply peaked to broad. The lead spectral densities are defined analogously as

[TABLE]

with both flat and Lorentzian forms considered below. Combining Eqs. (II.1), (II.1) and (II.1) we arrive at the weak-coupling BM master equation

[TABLE]

which on imposing the infinite bias limit becomes equivalent (for our spectral density) to the weak-coupling approaches considered for example in Refs. Stones and Olaya-Castro (2017); Kießlich et al. (2007). Our concern is primarily with finite bias, however.

II.2 Implementing the RC mapping

Having looked at the weakly-coupled master equation in the previous section we now consider the RC formalism and how it can be used to treat the same model though with strong electron-phonon coupling. To implement the RC mapping we apply a unitary normal mode transformation to our Hamiltonian, $H$, that incorporates a collective environmental degree of freedom (the RC) into the system (see Fig. 2), as previously applied to phonon coupling in other contexts Iles-Smith et al. (2014, 2016); Strasberg et al. (2016); Garg et al. (1985); Newman et al. (2017); Maguire et al. (2018). The strong coupling between our original system, S, and the RC results in a repartitioning of the original system and environment. As such we attain an augmented system, S*′, which includes both the original system and the RC, along with an explicit interaction between the two (denoted by $\kappa$). We treat the residual environment, E′*, as being weakly coupled to the RC and incorporate its influence using standard BM methods.

The RC approach has been successfully benchmarked against the numerically-exact hierarchical equations of motion over a range of parameters, including both broad and peaked spectral densities Iles-Smith et al. (2014, 2016) and situations in which system-environment correlations and non-Markovian effects are pronounced. For low temperatures, the residual environment correlation time may become extended such that the BM approximation on E*′* no longer applies. We do not consider such temperature regimes here, however.

On applying the RC mapping as described to $H$ of Eq. (1) we find

[TABLE]

Here, the augmented system Hamiltonian becomes

[TABLE]

where the first line contains terms associated with the original system, and we now also have its coupling to the RC as well as the self energy of the latter in the second line. The RC (collective coordinate) is defined directly from the electron-phonon interaction term in Eq. (5) such that

[TABLE]

where the RC creation ($a^{\dagger}$) and annihilation ($a$) operators satisfy $[a,a^{\dagger}]=1$, $\kappa^{2}=\sum_{q}g_{q}^{2}$ and the RC frequency is $\Omega=\omega_{0}$ Iles-Smith et al. (2014, 2016). The residual phonon environment E*′* is defined in terms of bosonic creation and annihilation operators $b_{k}^{\dagger}$ and $b_{k}$, respectively, for modes of frequency $\nu_{k}$,

[TABLE]

and is coupled with strengths $f_{k}$ directly to the RC by Iles-Smith et al. (2014, 2016); Hartmann et al. (2000)

[TABLE]

The RC mapping has no impact on the lead operators and so Eqs. (4) and (6) remain unchanged. We thus denote the total mapped Hamiltonian by

[TABLE]

Having outlined the mapping at the Hamiltonian level, we now consider the dynamical treatment of phonon interactions within the RC framework, before moving on to examine the lead couplings to our augmented system. The RC master equation has previously been derived in several different contexts Iles-Smith et al. (2014, 2016); Strasberg et al. (2016); Newman et al. (2017) and gives a Liouvillian of the form

[TABLE]

Here $\rho_{\rm S^{\prime}}(t)$ is the reduced density operator of the augmented system (comprising both S and RC), $A_{\rm ph}=(a^{\dagger}+a)$, and

[TABLE]

with $H_{\rm S^{\prime}}|\psi_{j}\rangle=\psi_{j}|\psi_{j}\rangle$, $A_{\mathrm{ph}_{jk}}=\langle\psi_{j}|A_{\rm ph}|\psi_{k}\rangle$, $\eta_{jk}=\psi_{j}-\psi_{k}$ and $J_{\rm RC}=(1/2\pi)\omega e^{-\omega/\Lambda}$. In practice, we represent the augmented system within a basis consisting of Fock states for the RC mode and the double quantum dot states $\{|0\rangle,|L\rangle,|R\rangle\}$. We take fourteen Fock states in all our calculations, which ensures numerical convergence. For the sake of completeness, we include a derivation of the RC Liouvillian within Appendix B.

We now turn our attention to the lead couplings and derive their contributions when coupled to our augmented system. It is common to simplify the model by taking the infinite (source-drain) lead bias limit, $\mu_{\rm L}=-\mu_{\rm R}\rightarrow\infty$, see e.g. Refs. Stones and Olaya-Castro (2017); Kießlich et al. (2007); Harbola et al. (2006b). In this regime, provided the lead spectral densities are flat, it is accurate to treat them by simply adding the standard Liouvillian terms to the master equation [see Eqs. (A) and (A)] unmodified due to the (potentially strong) phonon interactions. However, we wish to extend our considerations to finite bias regimes, in which case we shall see that including non-additive contributions becomes crucial. In order to do this we must derive the lead couplings with respect to our augmented system. In the following we examine the left lead contributions explicitly; the treatment of the right lead is analogous and so shall be stated afterwards. Setting $A^{\prime}_{1}=\ket{0}\bra{L}$ and $A^{\prime}_{2}=\ket{L}\bra{0}$ allows us to write the left lead interaction Hamiltonian, Eq. (6), in the form

[TABLE]

Here we have moved into the interaction picture with respect to the mapped system Hamiltonian $H_{\rm S^{\prime}}$, i.e. $A^{\prime}_{1(2)}(t)=e^{iH_{\rm S^{\prime}}t}A_{1(2)}^{\prime}e^{-iH_{\rm S^{\prime}}t}$, with $B_{1(2)}(t)=e^{iH_{\rm E_{L}}t}B_{1(2)}e^{-iH_{\rm E_{L}}t}$ as before. Now, following the standard BM approach Schaller (2014) within the mapped representation, we arrive at a Schrödinger picture Liouvillian of the form

[TABLE]

where $\rho_{\rm S^{\prime}}(t)$ is again the reduced density operator of the mapped (augmented) system comprising of both the DQD and the RC. As with the phonon case, we express the system operators within the eigenbasis of $H_{\rm S^{\prime}}$:

[TABLE]

Here, $\eta_{jk}$ and $\ket{\psi_{j}}$ are defined as before, with $A^{\prime}_{1_{j,k}}=\bra{\psi_{j}}A^{\prime}_{1}\ket{\psi_{k}}$ and $A^{\prime}_{2_{j,k}}=\bra{\psi_{j}}A^{\prime}_{2}\ket{\psi_{k}}$. We also consider the continuum limit for the lead coupling, taking $\Gamma_{\rm L}(\omega)=2\pi\sum_{k_{\rm L}}|t_{k_{\rm L}}|^{2}\delta(\omega-\epsilon_{k_{\rm L}})=2\pi J_{\mathrm{L}}(\omega)$ to give

[TABLE]

To simplify this expression we now consider a flat lead spectral density, such that $\Gamma_{\rm L}(-\eta_{jk})=\Gamma_{\rm L}(\eta_{jk})=\Gamma_{\rm L}$. We also define

[TABLE]

which allows the Liouvillian to be written in the more compact form

[TABLE]

The right lead Liouvillian can be given analogously as

[TABLE]

where $A^{\prime}_{3}=\ket{0}\bra{R}$, $A^{\prime}_{4}=\ket{R}\bra{0}$, and we define

[TABLE]

with $A^{\prime}_{3_{j,k}}=\bra{\psi_{j}}A^{\prime}_{3}\ket{\psi_{k}}$ and $A^{\prime}_{4_{j,k}}=\bra{\psi_{j}}A^{\prime}_{4}\ket{\psi_{k}}$. Finally we combine Eqs. (II.2), (II.2) and (II.2) to arrive at the full master equation

[TABLE]

Though this is additive at the level of the augmented system S*′* due to the BM approximations made for both the lead couplings and the residual phonon bath, it encodes non-additive contributions at the level of the original system $S$ due to the explicit presence of coupling to the RC and the resultant changes in the lead Liouvillians. On the other hand, a truly additive master equation is obtained following the RC mapping if we simply append the standard weak-coupling lead Liouvillians, Eqs. (II.1) and (II.1), to the phonon Liouvillian of Eq. (II.2):

[TABLE]

where $\openone_{\rm RC}$ is the identity operator in the RC subspace. This procedure implements the RC mapping for the phonons, but the lead dissipators are now unaware of the augmented system and so are treated phenomenologically. We shall see that this results in inaccuracies for flat leads when we consider finite bias regimes, and even at infinite bias for Lorentzian lead spectral densities.

II.3 Full counting statistics

Electron counting statistics provides a means to track the electron flow between our system and a reservoir of our choosing. In this work we consider the transfer of electrons between the right lead and the right dot, though we note that within the steady-state it would be equivalent to track the flow between the left lead and the left dot. As mentioned previously, the fact that the RC treatment results in a Markovian master equation in the mapped representation for our originally non-Markovian problem implies that the counting statistics formalism can be applied in the standard manner. We thus consider single electron transitions at each stage of the evolution, and as such we can take the description of a generic reduced density operator $\dot{\rho}(t)=\mathcal{L}\rho(t)$ with appropriate Liouvillian superoperator $\mathcal{L}$, and split it into three separate terms

[TABLE]

Here we explicitly consider the number of electrons present on the right lead during each transition, as denoted by the superscripts present on $\rho(t)$. The first term includes the cases in which no electron transfers between the right dot and the right lead occur, thus the number of electrons present on the lead, $n$, is unaffected. The remaining two terms consider the cases of electron movement onto and off the right lead, such that $\mathcal{I}^{+}$ is the rate at which the right lead gains electrons ($n-1\rightarrow n$) and $\mathcal{I}^{-}$ is the rate at which it loses them ($n+1\rightarrow n$).

While it is possible to compute the full probability distribution $P_{n}(t)=\mathrm{Tr}(\rho^{(n)}(t))$ of the electron tunneling events, it is generally computationally expensive and time consuming. It is often more useful to consider the cumulant generating function $\mathcal{C}(\Psi,t)$, which is defined as

[TABLE]

We can then calculate the cumulants of our charge distribution using $\mathcal{C}(\Psi,t)$. We focus specifically on the first and second order cumulants as these are proportional to the mean particle current and zero-frequency noise when taking time derivatives:

[TABLE]

Following a rigorous derivation discussed in Flindt et al. (2010) we are able to define current and noise in an alternative way:

[TABLE]

where $\mathcal{I}=\mathcal{I}^{+}-\mathcal{I}^{-}$, $\mathcal{J}=\mathcal{I}^{+}+\mathcal{I}^{-}$, $|0\rangle\rangle$ is the right eigenvector of our Liouvillian corresponding to the zero eigenvalue, and $\langle\langle\tilde{0}|$ is the identity operator (the left eigenvector of the zero eigenvalue). We have also introduced the pseudoinverse of our Liouvillian, $\mathcal{R}=\mathcal{QL}^{-1}\mathcal{Q}$, with projectors $\mathcal{P}=|0\rangle\rangle\langle\langle\tilde{0}|$ and $\mathcal{Q}=\openone-\mathcal{P}$. The pseudoinverse is well defined as the inversion is performed in the subspace spanned by $\mathcal{Q}$. For more details see Flindt et al. (2005, 2010).

We shall also examine the Fano factor of our system, $F$, defined as the ratio between the noise and the current, i.e. a comparison between the variance of the data and its mean:

[TABLE]

The Fano factor can be used to signal deviations of our current fluctuations from Poissonian statistics. For $F>1$ we have super-Poissonian behaviour, indicating that the transport channel is disordered. This means that there will be times whereby few electrons tunnel into the right lead, and then at other times there will be a rapid succession of electron tunneling events. For $F<1$ we have sub-Poissonian behaviour. In this case the tunneling of electrons between the dots and the right lead is very uniform, and as such electrons appear to tunnel at evenly spaced intervals. As a final point we note the changes to Eqs. (50) and (51) when we are in the infinite bias regime, in which case $\mathcal{I}=\mathcal{J}=\mathcal{I}^{+}$ so we have

[TABLE]

III Results and Discussion

III.1 Impact of phonon coupling

We begin our results discussion by considering the impact that the phonon coupling strength has on the overall current flowing into the right lead, when both leads are taken to be flat, ${\mathcal{J}}_{\rm L/R}(\omega)=\Gamma_{\rm L/R}/2\pi$. We vary the coupling strength by altering the phonon reorganisation energy $\lambda$. In Fig. 3(a) we consider the infinite bias limit, whereby we take the chemical potential difference of our leads to be so large that they become an effective source (left lead) and drain (right lead). As such we insist on unidirectional flow of electrons, preventing any movement of electrons onto the left lead or off the right lead. From this figure we see that for the infinite bias regime there are no non-additive effects, and as such the additive (i.e. phenomenological lead dissipator) RC treatment of Eq. (II.2) and the full non-additive master equation [Eq. (II.2)] provide complete agreement for any value of $\lambda$. The BM approach displays qualitatively different behaviour at strong reorganisation energies where the weak phonon coupling theory breaks down. As we increase $\lambda$ from zero we see that all three approaches predict a growth in current. In addition to the primary transport mechanism (characterised by the inter-dot tunneling $T$), for larger phonon coupling the system is able to make better use of the secondary phonon-assisted channel, and as such the current flow increases. For values of the phonon reorganisation energy that are the same order of magnitude as the dot energy levels we find that the current is enhanced compared to the BM predictions, indicating that only the strong coupling techniques are able to fully utilise both transport channels. For even larger values of $\lambda$ the current flow within the strong coupling theories then drops rapidly. In this regime the substantial RC displacements lead to a strong suppression of vibrational overlaps, which renormalises the effective inter-dot tunneling strength towards zero (not captured within the BM theory). This in turn inhibits electronic transitions and results in localisation of electrons on the left dot, such that no current can flow.

Moving away from the infinite bias regime we consider cases whereby we relax the stringent tunneling constraints that it imposes on our system. In Fig. 3(b) we examine a finite bias value where $\mu_{\rm L}=-\mu_{\rm R}=250$ meV. Here, the left lead chemical potential is still larger than the DQD energies, encouraging transport from left to right, however the previous source-drain picture no longer applies. Though increasing the phonon coupling from zero again results in an increased current flow, the additive treatment no longer appears to experience full localisation at very large reorganisation energies. In contrast, the non-additive approach again shows complete localisation for very large values of $\lambda$, however, unlike the non-additive case it no longer predicts any improvement in current over the weak-coupling BM master equation. It is important to note that within the non-additive treatment, the leads become sensitive to DQD energies that have been shifted by the inclusion of the RC, which is not the case for the additive approach. Here, this shift acts to increase backflow of electrons in comparison to the additive theory, and thus removes the region of improved current flow seen at infinite bias, or indeed at lower biases as in Fig. 3(c).

When considering Fig. 3(c) we note that the DQD eigenstate $\ket{\Upsilon_{1}}$ is now higher in energy than the left lead chemical potential and as such it falls outside the bias window. This results in tunnelling through only a single level, $\ket{\Upsilon_{2}}$, and thus an overall reduction in current flow. For the BM approach there is now no longer any improvement in current as we increase the phonon coupling, whereas the two strong coupling approaches demonstrate a slight increase for large reorganisation energies, with the non-additive treatment predicting this to occur at lower reorganisation energies than the additive theory. We also see that, as in Fig. 3(b), the additive approach does not display the same level of electron localisation as the non-additive theory at very large reorganisation energies, and so does not experience as great a suppression of current.

We attribute the differences in predictions between the non-additive and additive theories at large reorganisation energies to a breakdown of the additive approximation. This can be seen explicitly in the case of zero bias, as shown in Fig. 3(d). Here, we would expect no current flow, as both leads are held at equal chemical potentials and all three baths (leads and phonons) have the same temperature. As such, there should be no energy gradients for our electrons to follow. This is correctly depicted by the non-additive RC mapping and even the BM approach, though the latter does not more generally apply when the reorganisation energy becomes significant. However, the additive RC case displays markedly incorrect behaviour, predicting a current flow of over $1.5$ nA for large values of the phonon coupling. This results from the leads being insensitive to the strong phonon coupling within the additive treatment, which thus predicts a non-equilibrium steady state even at zero bias. Specifically, the true equilibrium state at strong-coupling should be defined with respect to the augmented system of the double dot plus RC. This is accurately captured within the full non-additive treatment, but within the additive approximation the leads couple directly to the dots without reference to the RC. They thus cannot act to equilibrate the system within the correct, augmented basis, which results here in the incorrect generation of a non-equilibrium steady state even in the absence of chemical potential and temperature gradients. The associated spurious current may be positive or negative depending on the energetics of the dots and the phonon coupling, which varies as the reorganisation energy is increased.

III.2 Impact of voltage bias

Having considered the impact that phonon coupling has on the current flow we now turn our attention to the lead chemical potential bias. We have already demonstrated that non-additive contributions can result in significant differences when we consider finite bias regimes, and we aim in this section to consolidate that knowledge and outline further interesting behaviour that can be missed with a phenomenological additive treatment.

We start, however, with Fig. 4(a), in which we consider weak phonon coupling (with both the leads and the phonon bath held at room temperature) and so expect to see complete agreement between all of our approaches. Nevertheless, this figure allows us to highlight a number of aspects of the physics of our DQD system. Firstly, we note that even when the bias is negative (encouraging right to left transitions) we obtain a lower current flow than in the forward tunneling case. This is a by product of the energy levels present within our DQD system. Although the coupling between each dot and their respective lead is the same, the difference in energies of our dots results in an asymmetry within the system.

Another key feature of the current-voltage graphs is their usefulness in determining at what bias values our energy levels $\ket{\Upsilon_{1}}$ and $\ket{\Upsilon_{2}}$ enter into the bias window. Considering positive bias, from Fig. 4(a) we are able to determine two plateau-like regions in the current. In the first, we have achieved the maximum current possible for tunneling through the $\ket{\Upsilon_{2}}$ state. For the parameters considered, this state is held at an energy of roughly $-10$ meV and so is always within the bias window, except around $V=0$. On the other hand, $\ket{\Upsilon_{1}}$ has a much larger energy of $220$ meV and is therefore initially inaccessible. When the bias voltage is increased towards $500$ meV we see that the current begins to increase again, as now both system eigenstates lie within the bias window. Further increasing the bias leads to an eventual plateau in the current, as the system is now fully utilising its transport channels. Qualitatively similar behaviour is seen in Fig. 4(b) where the lead and phonon temperatures have been reduced, resulting in sharper current features.

Also shown in Figs. 4(a) and (b) are the corresponding current-voltage curves in the absence of phonon coupling. For positive bias, we see that even a weak phonon coupling strength results in a significant increase in current. Again, this can be attributed to the opening up of the secondary phonon-assisted transport mechanism. For the temperatures considered here phonon emission dominates over absorption, and thus phonons mediate transitions from the higher lying DQD eigenstate $\ket{\Upsilon_{1}}$ to the lower lying state $\ket{\Upsilon_{2}}$, resulting in a current enhancement. For the same reason, current is slightly suppressed at negative bias, whereby the phonons preferentially mediate transitions in the opposite direction to the overall flow of current.

These features are even more pronounced at strong phonon coupling, as depicted in Figs. 4(c) and (d), where we now also see substantial variations between the different theoretical approaches. Consistent with our previous discussion, the additive approach is inaccurate for small bias values, again giving a finite current even for vanishing bias voltage, though the additive and non-additive theories do tend towards each other in large bias regimes. Furthermore, we can now see significant discrepancies in predictions for the locations of plateaus within the current curves, particularly between the non-additive technique and the others, as well as a systematic underestimation within the weak phonon coupling approach of the current at large positive bias, and a corresponding overestimation of the current magnitude at large negative bias. The non-additive theory also displays a suppression of current for small bias voltages that is not reproduced by either the additive or BM methods. This is due to phonon-induced shifts in the electronic energy levels at strong phonon coupling that have the effect of keeping them out of the bias window for a larger voltage range, and is only correctly captured within the non-additive theory due to the sensitivity of the leads to the augmented (mapped) system. As in Franck-Condon blockade Koch and von Oppen (2005b), the underlying cause for these shifts is phonon displacement, here captured dynamically through the RC. Likewise, we find a series of plateaus in the non-additive case for positive bias at lower temperatures in Fig. 4(d), as consecutive states of the augmented system enter the bias window with increasing voltage difference Strasberg et al. (2016). In contrast, the leads within the additive and BM theories remain sensitive only to the original DQD eigenstates even at strong phonon coupling, and so exhibit qualitatively similar behaviour to the weak phonon coupling curves.

III.3 Noise and Fano factor

We now move on to look at the impact of non-additive effects on higher order cumulants. We shall focus here on the second order cumulant, which in the case of quantum transport problems is proportional to the zero-frequency noise, as well as its ratio to the first order cumulant, i.e. the Fano factor. Considering first the infinite bias case, in Fig. 5(a) we plot the second order cumulant, $\langle\langle I^{2}\rangle\rangle$, as a function of phonon coupling strength for the same set of parameters as the current in Fig. 3. We see that though the overall trends are similar here between $\langle\langle I^{2}\rangle\rangle$ and current, in the region in which the current begins to increase, $\lambda\approx 0.01-0.1$ meV, the second order cumulant displays a small decline, leading to more ordered transport as indicated by the decrease in Fano factor in Fig. 5(b). As $\lambda$ is increased further $\langle\langle I^{2}\rangle\rangle$ then follows the same trend as the current as shown by the Fano factor levelling off. Finally, at very large $\lambda$, the localisation of the electron onto one of the dots results in a sharp decrease in $\langle\langle I^{2}\rangle\rangle$. Nevertheless, for the largest reorganisation energies considered the current decreases quicker than $\langle\langle I^{2}\rangle\rangle$, resulting in a transition from a sub-Poissonian to a super-Poissonian Fano factor. As we are working within the Coulomb blockade limit this increased Fano factor does not imply electron bunching within our system, but rather indicates the random nature of the tunneling events Emary et al. (2012). Qualitatively similar behaviour is observed at finite bias in Figs. 5(c) and (d), though now we see the expected discrepancies between the additive and non-additive theories at large phonon reorganisation energies, with the non-additive Fano factor never becoming super-Poissonian.

III.4 Lorentzian Leads

So far we have only considered the simplest case of leads with flat spectral densities, using $\Gamma_{\rm L}$ and $\Gamma_{\rm R}$ to denote their couplings to the respective dots. We now generalise our treatment to Lorentzian lead spectral densities,

[TABLE]

with widths $\delta_{\mathrm{L/R}}$ and peak positions $v_{\mathrm{L/R}}$. This definition results in non-additive contributions being present even in the infinite bias case, as can be seen in Fig. 6(a) for the current and (b) for the Fano factor. For flat spectral densities at infinite bias the leads become completely insensitive to the eigenstructure of the system Hamiltonian, which results mathematically from the lead correlation functions $C_{12}(\tau)$ and $C_{43}(\tau)$ reducing to delta functions in time ($C_{21}(\tau)=0=C_{34}(\tau)$ in the same limit). This is not true for Lorentzian lead spectral densities. Since the leads are then sensitive to different systems in the additive and non-additive theories, the original S and augmented S*′*, respectively, the agreement between the two that was present at infinite bias for flat leads no longer holds. It is also worth noting from comparison of Figs. 3(a) and 6(a) that the overall current is decreased in the Lorentzian lead case (for the same $\Gamma_{L}$ and $\Gamma_{R}$) due to the reduced effective coupling strength away from $\omega=v_{\rm L/R}$. Consequently, phonon coupling effects tend to be slightly more pronounced for Lorentzian leads, given the same set of parameters.

IV Conclusions

We have employed a combination of the RC mapping and counting statistics techniques to probe electron tunneling through a molecular nanojunction (DQD) system with strong and non-Markovian coupling to phonons. By deriving the lead contributions to the dynamics within the RC mapped (augmented) system representation, we are able to incorporate non-additive phonon-lead effects into the electron counting statistics formalism, which arise here due to the non-negligible phonon interactions.

We find that for flat lead spectral densities, non-additive effects become pronounced for finite lead bias voltages, and that a phenomenological additive approach in which the leads and phonons are treated completely separately can even give rise to unphysical features, such as non-zero current for vanishing lead bias (with all environments held at the same temperature). For Lorentzian leads, non-additive effects manifest themselves even within the infinite bias regime, due to the non-trivial energy dependence of the lead couplings.

In future work it will be interesting to consider the impact of non-additive contributions to the thermoelectric properties of molecular nanojunctions, where we no longer restrict ourselves to regimes in which the phonons and leads are held at the same temperature. It would also be intriguing to extend the approach to modelling strong lead couplings, by generalising the RC formalism to fermionic environments Strasberg et al. (2018); Nazir and Schaller (2018). A further avenue of investigation would be to implement counting statistics on the phonon bath to monitor heat flow between the system and the vibrational environment.

Acknowledgements.

We thank Neill Lambert and Gernot Schaller for discussions. C.M. is supported by the EPSRC. A.N. is supported by the EPSRC, grant no. EP/N008154/1.

Appendix A Derivation of Weakly Coupled Leads

When considering the lead contributions we transform the relevant terms into the system’s energy eigenbasis. To do so we diagonalise the system Hamiltonian as given in Eq. (2). Following the diagonalisation the system Hamiltonian is

[TABLE]

where the $\Upsilon_{i}$’s are the eigenvalues of $H_{\rm S}$ and the $\ket{\Upsilon_{i}}$’s are the corresponding eigenvectors, given as

[TABLE]

Here we have defined $\Delta=\frac{\epsilon_{\rm L}-\epsilon_{\rm R}}{2}$ and $\alpha=\Delta-\sqrt{\Delta^{2}+T^{2}}$ for brevity. We also rearrange our eigenvector equations in order to define the original states in terms of the eigenbasis. Rearranging Eqs. $(\ref{E1})$ and $(\ref{E2})$ we find

[TABLE]

Having diagonalised the system Hamiltonian we now consider the left lead terms present in the interaction Hamiltonian and express them within the system eigenbasis. Considering Eq. (6) the interaction Hamiltonian is given as

[TABLE]

in terms of the system eigenbasis. We define $A_{1}=\alpha\ket{0}\bra{\Upsilon_{2}}-T\ket{0}\bra{\Upsilon_{1}}$ and $A_{2}=\alpha\ket{\Upsilon_{2}}\bra{0}-T\ket{\Upsilon_{1}}\bra{0}$ and find the time dependent form of the two operators in the interaction picture to be

[TABLE]

The environmental contributions are the same as the standard fermionic correlation functions $C_{12}(t)=\sum_{k_{\rm L}}|t_{k_{\rm L}}|^{2}e^{i\epsilon_{k_{\rm L}}t}f_{\rm L}(\epsilon_{k_{\rm L}})$ and $C_{21}(t)=\sum_{k_{\rm L}}|t_{k_{\rm L}}|^{2}e^{-i\epsilon_{k_{\rm L}}\tau}(1-f_{\rm L}(\epsilon_{k_{\rm L}}))$ Schaller (2014). The Liouvillian for considering only the left lead contributions is then given as

[TABLE]

After some straightforward algebra we obtain an explicit form for this Liouvillian in the energy eigenbasis, which we can express instead in the site basis making use of Eqs. (59) and (60). We then arrive at the following cumbersome expression for the left lead contribution for a flat spectral density

[TABLE]

We note that there are no terms that would indicate movement from the ground state to the right dot state (i.e. terms of the form $\ket{R}\bra{0}\rho_{\rm S}(t)\ket{0}\bra{R}$) or vice versa, as expected in our DQD set up. We do see, however, that in the finite bias case we have coherence contributions between the left and right dots, generated by the inter-dot tunneling parameter, $T$, present in the system Hamiltonian. If we impose the infinite bias limit by taking the Fermi factor $f_{\rm L}\rightarrow 1$ these coherence terms disappear and we are left with

[TABLE]

which describes hopping of electrons from the left lead to the left dot, rather than to the DQD eigenstates.

The right lead contributions can be stated analogously:

[TABLE]

As before we can impose the infinite bias limit taking $f_{\rm R}\rightarrow 0$ and we find the Liouvillian reduces to

[TABLE]

Eqs. (A) and (A) are the same as those given, for example, in Ref. Stones and Olaya-Castro (2017) for the lead contributions within the infinite bias regime.

Appendix B Electron-Phonon Coupling in the RC Formalism

When considering molecular nanojunction models we want to be able to model the impact of strong vibrational coupling on the system. In order to do this we make use of the RC mapping to incorporate a collective mode of the vibrational environment within an augmented system. The procedure to do so is discussed here. We ignore the leads for simplicity within this appendix, but their inclusion either additively or non-additively is discussed in the main text.

We begin by considering the post-mapping interaction Hamiltonian of Eq. (24). Defining $A_{\mathrm{ph}}=(a^{\dagger}+a)$, $B_{\mathrm{ph}}=\sum_{k}f_{k}(b_{k}^{\dagger}+b_{k})$ and $\xi=\sum_{k}\frac{f_{k}^{2}}{\nu_{k}}$ we can write

[TABLE]

We now derive the associated Liouvillian, beginning with the Liouville-von Neumann equation

[TABLE]

for the full system-environment state $\tilde{\rho}(t)$ within the interaction picture with respect to $H_{\rm S^{\prime}}+H_{\rm E^{\prime}}$. Treating the residual bath with the BM approximations we arrive at the second-order perturbative form

[TABLE]

where $\tilde{\rho}_{\rm S^{\prime}}(t)$ is the reduced system-RC density operator within the interaction picture. Here we have defined the bath correlation functions $\Gamma^{\pm}(\tau)=\mathrm{Tr_{E}}\{(B_{\mathrm{ph}}(\tau)B_{\mathrm{ph}}\pm B_{\mathrm{ph}}(-\tau)B_{\mathrm{ph}})\rho_{\rm E^{\prime}}\}/2$, with $\rho_{E^{\prime}}$ a (stationary) thermal state of the residual environment. We can further simplify this Liouvillian by noting that $-i\xi[A_{\mathrm{ph}}^{2}(t),\rho_{\rm S^{\prime}}(t)]=-i\xi[A_{\mathrm{ph}}^{2}(t),\rho(0)]-\xi^{2}\int^{\infty}_{0}d\tau[A_{\mathrm{ph}}^{2}(t),[A_{\mathrm{ph}}^{2}(t-\tau),\rho_{\rm S^{\prime}}(t)]]$ allowing us to make a substitution that reduces the first two terms into a single commutator. Introducing the forms of the rates $\Gamma^{\pm}$ in the continuum limit,

[TABLE]

we are able to rewrite Eq. (B) as

[TABLE]

We take the residual bath spectral density to be $J_{\rm RC}(\omega)=\frac{1}{2\pi}\omega e^{-\omega/\Lambda}$ and consider the limit of infinite cutoff frequency ($\Lambda\rightarrow\infty$). We then obtain what appear to be divergent terms within our master equation, however it is possible to eliminate these terms by integrating by parts Iles-Smith et al. (2014, 2016). As such we can cancel out the divergences and upon moving back to the Schrödinger picture we have

[TABLE]

with $H_{\rm S^{\prime}}$ the mapped system-RC Hamiltonian as before. We now need to consider our mapped system operators, $A_{\mathrm{ph}}(-\tau)$. To proceed we express them within the augmented system eigenbasis as

[TABLE]

where the augmented system eigenequation $H_{\rm S^{\prime}}\ket{\psi_{j}}=\psi_{j}\ket{\psi_{j}}$ is solved numerically, and $A_{\mathrm{ph}_{j,k}}=\bra{\psi_{j}}A_{\mathrm{ph}}\ket{\psi_{k}}$. Within the interaction picture with respect to $H_{\rm S^{\prime}}$ we then have

[TABLE]

with $\eta_{jk}=\psi_{j}-\psi_{k}$ being the eigenvalue differences. Using these definitions we arrive at the final Liouvillian form

[TABLE]

Here we have neglected the imaginary Lamb shift terms and included the rates from Eq. (B) within the system operators:

[TABLE]