Photoinduced Anomalous Coulomb Blockade and the Role of Triplet States in Electron Transport through an Irradiated Molecular Transistor II: Effects of Electron-Phonon Coupling and Vibrational Relaxation

TL;DR

This paper extends a theoretical model to study how electron-phonon interactions and vibrational relaxation affect photoinduced anomalous Coulomb blockade in molecular transistors, revealing new photocurrent mechanisms due to these couplings.

Contribution

It introduces a generalized master equation approach using Redfield theory to analyze the effects of electron-phonon coupling and vibrational relaxation on Coulomb blockade phenomena.

Findings

Identification of a new photocurrent mechanism due to electron-phonon coupling.

Analytical formulas for photoinduced current in different vibrational regimes.

Predictions of observable photoinduced Coulomb blockade effects in experiments.

Abstract

We generalize our previous theory [Nano Lett. \textbf{18}, 5015-5023 (2018)] to investigate the influence of electron-phonon (e-p) coupling and vibrational relaxation on photoinduced \textit{anomalous Coulomb blockade}, which originates from the triplet states and the energy level alignment. We derive the master equation for an irradiated molecular transistor and obtain the relevant rates via the Redfield theory instead of the phenomenological Fermi golden rule approach. To explore the interplay between e-p coupling and vibrational relaxation, we analyze the charge stability diagrams and the current-voltage characteristics (both gate voltage and source-drain bias voltage) under different e-p coupling strengths in two extreme limits of vibrational relaxation (equilibrated and unequilibrated phonon regimes). From the perspective of energy level alignment, we choose four representative situations and derive the analytical formulas of the photoinduced current in the equilibrated regime. The analytical solution reveals a new type of photocurrent due to e-p coupling that does not require the perfect energy level alignment between charged states and triplet states. In general, our study indicates that photoinduced current and anomalous Coulomb blockade caused by the triplet states are supposed to be experimentally observed.

Tables

Table 1. Table 1: Parameters adopted in the computations. Details are discussed in main text.

Parameter		Value		Meaning
$E_{0,S_0}$		$0$ eV		Energy level of the singlet ground state $|\mathrm{\hspace{0.17em}0},S_0⟩$
$E_{0,S_1}$		$2.01$ eV		Energy level of the singlet excited state $|\mathrm{\hspace{0.17em}0},S_1⟩$
$E_{0,T_1^m}$		$1.00$ eV		Energy level of the degenerate triplet states $|\mathrm{\hspace{0.17em}0},T_1^m⟩$
$E_{-1,D_0^{\sigma }}$		$6.045$ eV		Energy level of the cation doublet states $|-1,D_0^{\sigma }⟩$, image charge correction included.
$E_{+1,D_0^{\sigma }}$		$-2.06$ eV		Energy level of the anion doublet states $|+1,D_0^{\sigma }⟩$, image charge correction included.
$\hslash {\omega }_{\mathrm{vib}}$		$0.2$, $0.02$ eV		Energy of the intramolecular vibrational mode
$\lambda$		$0$,$0.25$,$0.5$,$1.0$,$2.0$		Electron-Phonon coupling
${\mu }_0$		$-5.3$ eV		Fermi energy in the electrode at zero source-drain bias voltage
$\hslash \Gamma$		${10}^{-4}$ eV		${\Gamma }_S={\Gamma }_D=\Gamma$ in a symmetric molecular junction
$k_{B}T$		0.05$\hslash {\omega }_{\mathrm{vib}}$		Thermal energy at temperate T
$k_{0,T_1\leftarrow 0,S_1}$		${10}^6$ $s^{-1}$		Intersystem crossing rate
$k_{0,S_0\leftarrow 0,T_1}$		${10}^3$ $s^{-1}$		Phosphorescence rate
$k_{0,S_0\leftarrow 0,S_1}^{\mathrm{spon}}$		${10}^8$ $s^{-1}$		Spontaneous emission rate

Equations

$${\hat{H}}_{\mathrm{tot}}={\hat{H}}_0+\hat{V},$$

$${\hat{H}}_0={\hat{H}}_{\mathrm{mol}}+{\hat{H}}_{\mathrm{lead}}+{\hat{H}}_{\mathrm{th}}.$$

$$\hat{V}={\hat{H}}_{\mathrm{SOC}}+{\hat{H}}_{\mathrm{field}}+{\hat{H}}_{\mathrm{m}-\mathrm{th}}+{\hat{H}}_{\mathrm{m}-\mathrm{l}},$$

$$\displaystyle\hat{H}_{\mathrm{lead}}=\sum_{\alpha}\hat{H}_{\mathrm{lead,\alpha}}=\sum_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma}\varepsilon_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma}\hat{a}^{\dagger}_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma}\hat{a}_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma},$$

$$\displaystyle\hat{\rho}_{\mathrm{lead}}=\prod_{\alpha=S,D}\frac{\mbox{$\displaystyle\mathbf{e}^{-\beta\mbox{$\displaystyle\left(\hat{H}_{\mathrm{lead,\alpha}}-\mu_{\alpha}\hat{N}_{\mathrm{lead},\alpha}\right)$}}$}}{\mathrm{Tr}_{\mathrm{lead}}\{\mbox{$\displaystyle\mathbf{e}^{-\beta\mbox{$\displaystyle\left(\hat{H}_{\mathrm{lead,\alpha}}-\mu_{\alpha}\hat{N}_{\mathrm{lead},\alpha}\right)$}}$}\}},$$

$$\displaystyle\hat{H}_{\mathrm{th}}=\sum_{\alpha}\mbox{$\displaystyle\left(\frac{\hat{p}^{2}_{\alpha}}{2m_{\alpha}}+\frac{1}{2}m_{\alpha}\omega^{2}_{\alpha}\hat{x}^{2}_{\alpha}\right)$},$$

$$\displaystyle\hat{H}_{\mathrm{th}}=\sum_{\alpha}\mbox{$\displaystyle\left(\hat{b}^{\dagger}_{\alpha}\hat{b}_{\alpha}+\frac{1}{2}\right)$}\hbar\omega_{\alpha},$$

$$\displaystyle\hat{\rho}_{\mathrm{th}}=\mbox{$\displaystyle\mathbf{e}^{-\beta\hat{H}_{\mathrm{th}}}$}/\mathrm{Tr}_{\mathrm{th}}\{\mbox{$\displaystyle\mathbf{e}^{-\beta\hat{H}_{\mathrm{th}}}$}\}.$$

$${\hat{H}}_{\mathrm{mol}}={\hat{H}}_{\mathrm{m}-\mathrm{el}}+{\hat{H}}_{\mathrm{m}-\mathrm{vib}}.$$

$${\hat{H}}_{\mathrm{m}-\mathrm{el}}$$

$$\displaystyle\hat{H}_{\mathrm{m-vib}}=\sum_{|N,a}\biggr{[}\frac{\hat{p}^{2}}{2m_{\mathrm{vib}}}+\frac{1}{2}m_{\mathrm{vib}}\omega^{2}_{\mathrm{vib}}\mbox{$\displaystyle\left(\hat{q}-q_{N,a}\right)$}^{2}\biggr{]}|N,a\rangle\langle N,a|,$$

$${\hat{H}}_{\mathrm{m}-\mathrm{vib}}=$$

$$\displaystyle-\lambda_{N,a}\mbox{$\displaystyle\left(\hat{b}^{\dagger}+\hat{b}\right)$}+\lambda_{N,a}^{2}\biggr{]}\mbox{$\displaystyle|\,N,a\,\rangle$}\mbox{$\displaystyle\langle\,N,a\,|$},$$

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

$$\displaystyle+V^{\mathrm{SOC}}_{0T^{m}_{1},0S_{1}}\mbox{$\displaystyle|\,0,T^{m}_{1}\,\rangle$}\mbox{$\displaystyle\langle\,0,S_{1}\,|$}\biggr{)},$$

$$H_{\mathrm{field}}$$

$$\displaystyle\mbox{$\displaystyle\textbf{E}$}\mbox{$\displaystyle\left(t\right)$}=\mbox{$\displaystyle\textbf{E}$}_{0}\mbox{$\displaystyle\left(\mbox{$\displaystyle\mathbf{e}^{i\omega t}$}+\mbox{$\displaystyle\mathbf{e}^{-i\omega t}$}\right)$}$$

$$\displaystyle\hat{H}_{\mathrm{m-th}}=\sum_{N,a}\sum_{\alpha}c_{\alpha}\mbox{$\displaystyle\left(\hat{q}-q_{N,a}\right)$}\hat{q}_{\alpha}|N,a\rangle\langle N,a|,$$

$${\hat{H}}_{\mathrm{m}-\mathrm{th}}={\hat{S}}_{\mathrm{m}-\mathrm{th}}{\hat{B}}_{\mathrm{m}-\mathrm{th}},$$

$$\displaystyle\hat{S}_{\mathrm{m-th}}=\sum_{N,a}\mbox{$\displaystyle\left(\hat{q}-q_{N,a}\right)$}|N,a\rangle\langle N,a|$$

$${\hat{B}}_{\mathrm{m}-\mathrm{th}}=\sum _{\alpha }{c}_{\alpha }{\hat{q}}_{\alpha }.$$

$${\hat{H}}_{\mathrm{m}-\mathrm{l}}$$

$$+h.c.),$$

$$\displaystyle\hat{H}_{\mathrm{m-l}}=\sum_{\alpha}\mbox{$\displaystyle\left(\hat{B}^{+}_{\alpha}\hat{S}^{-}_{\alpha}+\hat{S}^{+}_{\alpha}\hat{B}^{-}_{\alpha}\right)$}\,,$$

$$\displaystyle\hat{B}^{+}_{\alpha}=\mbox{$\displaystyle\left(\hat{B}^{-}_{\alpha}\right)$}^{\dagger}=\sum_{\mbox{$\displaystyle\textbf{k}$}\sigma}M_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma}\hat{a}_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma},$$

$$\displaystyle\hat{S}^{+}_{\alpha}=\mbox{$\displaystyle\left(\hat{S}^{-}_{\alpha}\right)$}^{\dagger}=\sum_{N,a,b}T^{\alpha}_{Na,N-1b}\mbox{$\displaystyle|\,N,a\,\rangle$}\mbox{$\displaystyle\langle\,N-1,b\,|$},$$

$$\displaystyle V_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma,NaN-1b}=T^{\alpha}_{NaN-1b}M_{\alpha\mbox{$\displaystyle\textbf{k}$}\sigma}\,.$$

$$\displaystyle\mbox{$\displaystyle|\,N,a,\nu\,\rangle$}=\mbox{$\displaystyle|\,N,a\,\rangle$}\otimes\mbox{$\displaystyle|\,\nu\mbox{$\displaystyle\left(\lambda_{N,a}\right)$}\,\rangle$},$$

$${\hat{H}}_{\mathrm{mol}}^{\prime }=$$

$$\displaystyle\sum_{N,a,\nu}\mbox{$\displaystyle\left(E_{N,a}-NeV_{G}+(\nu+\frac{1}{2})\hbar\omega_{\mathrm{vib}}\right)$}\mbox{$\displaystyle|\,N,a,\nu\,\rangle$}\mbox{$\displaystyle\langle\,N,a,\nu\,|$}\,,$$

$$\displaystyle\hat{U}=\mbox{$\displaystyle\mathbf{e}^{\sum_{N,a}\lambda_{N,a}\mbox{$\displaystyle|\,N,a\,\rangle$}\mbox{$\displaystyle\langle\,N,a\,|$}\mbox{$\displaystyle\left(\hat{b}^{\dagger}-\hat{b}\right)$}}$}\,.$$

$${\hat{H}}_{\mathrm{m}-\mathrm{th}}^{\prime }={\hat{S}}_{\mathrm{m}-\mathrm{th}}^{\prime }{\hat{B}}_{\mathrm{m}-\mathrm{th}}^{\prime }$$

Full text

††preprint: AIP/123-QED

I Introduction

Molecular electronics attracts broad interest in the field of chemical physics because it bridges fundamental sciences and frontier nanotechnology.Aviram and Ratner (1974); Nitzan and Ratner (2003); Joachim and Ratner (2005); Venkataraman et al. (2006); Ratner (2013); Huang et al. (2015); Hsu et al. (2017); Xin et al. (2019); Liu et al. (2019) During the past two decades, electron transport through a variety of molecules has been extensively investigated, and a lot of intriguing physical phenomena has been experimentally observed. Reed et al. (1997); Cui et al. (2001); Park et al. (2002); Xu and Tao (2003); Quek et al. (2009); Parks et al. (2010); Vazquez et al. (2012); Burzuri et al. (2014); Yoshida, Shibata, and Hirakawa (2015); Garner et al. (2018); Zhou et al. (2018); Lumbroso et al. (2018); Huang et al. (2019). The experimental breakthrough together with theoretical advances has made molecular electronics thriving and robust. However, for a practical device, the capability of tuning electric current by external fields is crucial and necessary. Thus, the manipulation of electron transport, e.g., by means of a back gateHeersche et al. (2006); Steele et al. (2009); Osorio et al. (2010); Perrin et al. (2013), electrochemical gateXiao et al. (2005); Huang et al. (2015), or lightJia et al. (2016); Fung et al. (2017); Park et al. (2002), is a primary focus for the whole field of molecular electronics.

Electric current through a molecular junction can be controlled by an optical field via either chemical or physical processes. Chemically, the high- and low-conductance states based on different molecular conformations have been able to be switched by a laser field at specific wavelengthsKatsonis et al. (2006); Dulić et al. (2003); Roldan et al. (2013); Jia et al. (2016). Physically, photoelectric current can be generated via the mechanisms of photon-assisted tunnelingSchiffrin et al. (2012); Paasch-Colberg et al. (2014); Yoshida, Shibata, and Hirakawa (2015); Cocker et al. (2016); Rybka et al. (2016); Yoshioka et al. (2016); Fung et al. (2017) or electronic excitationYasutomi et al. (2004); Banerjee et al. (2010); Battacharyya et al. (2011); Fainberg and Seideman (2012); Kornbluth, Nitzan, and Seideman (2013); Zhou et al. (2018); Morteza Najarian, Bayat, and McCreery (2018); Najarian and McCreery (2019) in a molecular junction. From the theoretical aspects, the Floquet-based methods have successfully demonstrated several fascinating phenomena, including photon-assisted tunnelingTien and Gordon (1963); Keller et al. (2002); Urdaneta et al. (2004, 2007); Tikhonov, Coalson, and Dahnovsky (2002); Hsu and Rabitz (2012); Hsu, Xie, and Rabitz (2014); Hsu and Rabitz (2015a), coherent destruction of tunnelingGrossmann et al. (1991); Camalet et al. (2003), coherent revival of tunnelingHsu and Rabitz (2015b), and quantum ratchet effectLehmann et al. (2002). For a single-photon process, kinetic approaches, e.g., quantum master equation approachesMuralidharan, Ghosh, and Datta (2006); Segal (2006); May and Kühn (2008a); Fainberg, Jouravlev, and Nitzan (2007); Zelinskyy and May (2011); Fainberg and Seideman (2012) or non-equilibrium green’s function (NEGF) approachesGalperin and Nitzan (2005); Galperin, Nitzan, and Ratner (2006); Galperin and Nitzan (2006); Galperin, Nitzan, and Ratner (2008); Chen, Ochoa, and Galperin (2017); Galperin (2017), offer a clear description of electron current induced by electronic excitation and optical excitation. However, in these theoretical studies, molecular triplet states are not considered in the modelling because of the slow rate of an intersystem crossing process. In fact, our previous studyFu et al. (2018) has clearly shown that the presence of the triplet states is significant for the light-driven transport properties regardless of whether the rate of the intersystem crossing is fast or not.

The triplet states play a crucial role in light-driven electron transport. In the previous studyFu et al. (2018), we introduce the concept of renormalized state energy and find that the low bias photoinduced current can be maximally achieved when "the renormalized state energy of the charged states is aligned with the energy of the triplet states", in which the triplet dark states indeed dominate the transport properties. Furthermore, our study shows that an optical field enables a field-off Coulomb diamond to be decomposed into three smaller diamonds, and the presence of the triplet states is revealed by the central diamond in the field-on conductance spectra. The concepts of renormalized state energy, energy-level alignment and photoinduced anomalous Coulomb blockade are general and allow us to explore the importance of dark states in light-driven quantum transport. However, the previous study does not consider the effects of electron-phonon (e-p) coupling and vibrational relaxation, and these two effects may interrupt the experimental observation of anomalous Coulomb blockade or photoinduced current attributed to the energy-level alignment of the triplet states and the charged states. In order to eliminate these concerns, a general theory of light-driven electron transport involving the triplet states, e-p coupling and vibrational relaxation is required.

Electron-vibration interactions exhibit rich physical phenomena in nanoscale electron transport and inspire extensive theoretical studies.Segal et al. (2000); Mitra, Aleiner, and Millis (2004); Wegewijs and Nowack (2005); Galperin, Ratner, and Nitzan (2005); Ryndyk, Hartung, and Cuniberti (2006); Hsu, Tsai, and Jin (2010a); Galperin, Ratner, and Nitzan (2007); Koch, Raikh, and von Oppen (2005); Koch and von Oppen (2005); Koch, von Oppen, and Andreev (2006); May and Kühn (2008a, b, c); Wang and May (2010a, b); Simine and Segal (2014); Agarwalla, Jiang, and Segal (2015) Among these studies, we would like to emphasize the pioneering works by KochKoch et al. (2004); Koch, Raikh, and von Oppen (2005); Koch and von Oppen (2005); Koch, von Oppen, and Andreev (2006) and MayMay and Kühn (2008a, b, c); Wang and May (2010a, 2011a, b, 2011b). Koch et. al. first predicted Franck-Condon blockadeKoch and von Oppen (2005) using a Pauli master equation approach based on the Fermi golden rule rates, where the effect of the vibrational relaxation is introduced phenomenologically via the single-mode relaxation time approximation. May et. al. discussed the photoinduced removal of the Franck-Condon blockade using a generalized rate equation approach.May and Kühn (2008a, b, c) The generalized rate equation approach is derived based on the electron-vibrational states using the projection operator method which includes the 2nd order rates of charge transfer transitions, vibrational relaxations, optical transitions and the molecular deexcitation due to the electrode. Based on the studies above, we present a Pauli master equation approach involving relevant transitions together with a pedagogical derivation without using the projection operator method.

This article is divided into five parts. In section II, we introduce our model (including a molecule, two leads, thermal bath, spin-orbit coupling (SOC), and light-matter interactions) and the corresponding Hamiltonians. To adequately treat the vibronic coupling, we apply the polaron transformation and analyze the model Hamiltonians in the polaron frame. In section III, in the framework of the Redfield theory, we derive the rates of charge transfer transition, vibration relaxation and optical transition together with phenomenologically including the rate of the intersystem crossing. Finally, we arrive at the master equation and the current formula in our model. In section IV, we discuss the rationality of the parameters in the work. In section V, first, we define the equilibrated and unequilibrated phonon regimes. Based on the two regimes, we explore the charge stability diagrams and the current-voltage curves of an irradiated molecular junction. Besides, we study the influence of vibronic coupling and bath relaxation on photoinduced Coulomb diamonds and anomalous Coulomb blockade. Moreover, according to the energy level alignment, we analyze the origin of the maximal photoinduced current and the unexpected photoresistive behavior.

II Model and Theory

II.1 Hamiltonian

As illustrated in Figure 1, we study the light-driven electron transport through a three-terminal molecular device, where a molecule is connected to a source (S) electrode and a drain (D) electrode under a gate electric field. We particularly focus on electric current due to the mechanism of incoherent sequential tunneling in the weak coupling limit, which can be reached by using a molecule with saturated alkyl linkersPark et al. (2002). The total Hamiltonian of the device is written as

[TABLE]

where

[TABLE]

$\hat{H}_{\mathrm{mol}}$ represents the isolated molecular system, $\hat{H}_{\mathrm{lead}}$ refers to the isolated leads (electrodes), $\hat{H}_{\mathrm{th}}$ describes the thermal bath, and the coupling Hamiltonian is

[TABLE]

which consists of the SOC $\hat{H}_{\mathrm{SOC}}$, the coupling to the radiation field $\hat{H}_{\mathrm{field}}$, the coupling to the thermal bath $\hat{H}_{\mathrm{m-th}}$ and the molecule-lead coupling between the molecule and the electrodes $\hat{H}_{\mathrm{m-l}}$.

The two electrodes including a source lead and a drain lead are described by a non-interacting electron gas model as

[TABLE]

where \hat{a}^{\dagger}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma} (\hat{a}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma}) creates (annihilates) an electron with momentum k, spin $\sigma$ and energy \varepsilon_{\alpha\mbox{\displaystyle\textbf{k}}\sigma} in the electrode $\alpha$ (S for source lead, and D for source lead). The electronic relaxation in the electrodes is assumed to be fast compared to all the dynamical processes on the molecular system. Therefore, the electrons in the electrodes are always at the thermal equilibrium. The thermal equilibrium density operator of the isolated electrodes reads

[TABLE]

where $\beta=1/k_{B}T$ with the Boltzmann constant $k_{B}$, \hat{N}_{\mathrm{lead,\alpha}}=\sum_{\mbox{\displaystyle\textbf{k}}\sigma}\hat{a}^{\dagger}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma}\hat{a}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma} is the number operator of electrode $\alpha$, and $\mu_{\alpha}$ denotes the chemical potential of electrons in the electrode $\alpha$. A source-drain bias $V_{\mathrm{SD}}$ applied symmetrically across the molecular junction is then taken into account by setting $\mu_{\mathrm{S}}=\mu_{0}+V_{\mathrm{SD}}/2$ and $\mu_{\mathrm{D}}=\mu_{0}-V_{\mathrm{SD}}/2$, where $\mu_{0}$ is the chemical potential of the two electrodes at zero bias.

The inactive intramolecular vibrations with respect to the electronic transitions, together with the environment, are treated as the thermal bath and then modelled as a set of harmonic oscillators in $\hat{H}_{\mathrm{th}}$,

[TABLE]

where $\hat{p}_{\alpha}$, $\hat{x}_{\alpha}$, $m_{\alpha}$ and $\omega_{\alpha}$ are respectively the momentum operator, the position operator, the mass and the vibrational frequency of the thermal bath mode $\alpha$. The thermal bath Hamiltonian is diagonalized by introducing \hat{q}_{\alpha}=\mbox{\displaystyle\left(\hat{b}^{\dagger}{\alpha}+\hat{b}{\alpha}\right)}\ell_{\alpha}/\sqrt{2} and \hat{p}_{\alpha}=i\mbox{\displaystyle\left(\hat{b}^{\dagger}{\alpha}-\hat{b}{\alpha}\right)}\hbar/(\sqrt{2}\ell_{\alpha}) with $\ell_{\alpha}=\sqrt{\hbar/(m_{{}_{\alpha}}\omega_{{}_{\alpha}})}$, which results in

[TABLE]

where $\hat{b}^{\dagger}_{\alpha}$ ($\hat{b}_{\alpha}$) creates (annihilates) a vibrational quantum of mode $\alpha$ with vibrational frequency $\omega_{\alpha}$. Similar to the electrons in the electrodes, the relaxation process of the phonon bath is also assumed to be faster than all relevant processes in the molecular system. Therefore, the thermal bath is always kept at its thermal equilibrium characterized by the density operator

[TABLE]

The molecular Hamiltonian consists of an electronic part and a vibronic part, i.e.

[TABLE]

The electronic Hamiltonian of the isolated molecule is expressed in terms of many-electron states that span the Fock space ${\cal F}$, which reads

[TABLE]

The $\displaystyle|\,N,a\,\rangle$ labels a many-electron state with $N$ referring to the number of excess electrons with respect to the neutral state of the molecule and $a$ labelling the electronic state for its spin multiplicity and energy level, which is consistent with the convention in the community of photophysics. The energy level of the state $\displaystyle|\,N,a\,\rangle$ that includes the energy shift caused by the image charge effect is denoted by $E_{N,a}$. The effect of the gate voltage $V_{G}$ is considered by $NeV_{\mathrm{G}}$. Although the Fock space is constructed by a complete set of many-electron states, we employ a truncated Fock space which includes a finite number of states that participate in the relevant electronic transitions. Our truncated Fock space is spanned by cation doublet charge states ($\displaystyle|\,-1,D^{\sigma=\pm 1/2}_{0}\,\rangle$, anion doublet and $\displaystyle|\,1,D^{\sigma=\pm 1/2}_{0}\,\rangle$), neutral singlet ground state ($\displaystyle|\,0,S_{0}\,\rangle$, neutral singlet excited state $\displaystyle|\,0,S_{1}\,\rangle$) and three neutral triplet states ($\displaystyle|\,0,T^{m=0,\pm 1}_{1}\,\rangle$), where the degeneracy of doublet and triplet states are considered.

We take a molecular vibration mode as the reaction coordinate and write the vibronic part of the molecular Hamiltonian as

[TABLE]

where $m_{\mathrm{vib}}$, $\omega_{\mathrm{vib}}$, $\hat{p}$ and $\hat{q}$ are respectively the effective mass, the frequency, the momentum operator and the position operator of the active vibrational motion, and $q_{N,a}$ is the equilibrium position of the potential energy surface (PES) associated with the electronic state $\displaystyle|\,N,a\,\rangle$. We introduce \hat{q}=\mbox{\displaystyle\left(\hat{b}^{\dagger}+\hat{b}\right)}\ell/\sqrt{2} and \hat{p}=i\mbox{\displaystyle\left(\hat{b}^{\dagger}-\hat{b}\right)}\hbar/(\sqrt{2}\ell) with $\ell=\sqrt{\hbar/(m_{\mathrm{vib}}\omega_{\mathrm{vib}})}$, and then transform the Hamiltonian $\hat{H}_{\mathrm{m-vib}}$ into

[TABLE]

where $\hat{b}^{\dagger}$ ($\hat{b}$) creates (annihilates) a vibronic quantum with frequency $\omega_{\mathrm{vib}}$, and the dimensionless e-p coupling parameter $\lambda_{N,a}$ is $q_{N,a}\sqrt{m_{\mathrm{vib}}\omega_{\mathrm{vib}}/2\hbar}$.

The coupling Hamiltonian $\hat{V}$ describes the couplings between the molecular system and the other external degrees of freedom, which are respectively responsible for the singlet-triplet transition, the radiation-induced transition, the vibrational relaxation and the charge transfer transition. The SOC (the origin of the singlet-triplet transitions) is represented as

[TABLE]

where the transition amplitudes $V^{\mathrm{SOC}}_{0S_{0},0T^{m}_{1}}$ and $V^{\mathrm{SOC}}_{0T^{m}_{1},0S_{1}}$ could be introduced by a relativistic correction to the kinetic terms. $\hat{H}_{\mathrm{SOC}}$ is responsible for the phenomenological rates of singlet-triplet transitions in the Pauli master equation. The explicit expression of the SOC and its originKöhler and Bässler (2009); Marian (2012); Powell (2015); Baryshnikov, Minaev, and Ågren (2017) are not discussed here.

The spin-allowed transitions between the neutral singlet states $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,0,S_{1}\,\rangle$, due to the coupling to a radiation field, are described by $\hat{H}_{\mathrm{field}}$ based on the semi-classical radiation theory and the dipole approximationSchatz and Ratner (1993) as

[TABLE]

where h.c. stands for Hermitian conjugate, \mbox{\displaystyle\textbf{d}}_{01} is the transition dipole moment between $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,0,S_{1}\,\rangle$, and E$\displaystyle\left(t\right)$ the electric field component of the radiation field. The time-dependent electric field takes the form of

[TABLE]

for a continuous wave laser source.

The vibrational relaxation process within a PES associated with the electronic state $\displaystyle|\,N,a\,\rangle$ originates from the coupling of the molecular vibrations to the thermal bath, i.e., $\hat{H}_{\mathrm{m-th}}$, which takes the following bilinear formCaldeira and Leggett (1983); Garg, Onuchic, and Ambegaokar (1985); Leggett et al. (1987); Jean, Friesner, and Fleming (1992); May and Kühn (2008a)

[TABLE]

where $c_{\alpha}$ characterizes the coupling the active molecular vibration to the thermal bath mode $\alpha$. For convenience of the derivation of the master equation, we write Eq. (16) in terms of a system oporator and a bath operator as

[TABLE]

where the system operator and the bath operator are respectively given by

[TABLE]

and

[TABLE]

Finally, we model the molecule-lead coupling in the following form,

[TABLE]

which leads to the exchange of electrons between the molecule and the two electrodes, and this coupling is the origin of the current through the molecular junction. The coupling parameter V^{*}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b} (V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b}) is a scattering amplitude which describes the tunneling of an electron from the molecule (electrode $\alpha$) to the electrode $\alpha$ (molecule) together with a transition of the molecule from $\displaystyle|\,N,a\,\rangle$ to $\displaystyle|\,N-1,b\,\rangle$ (from $\displaystyle|\,N-1,b\,\rangle$ to $\displaystyle|\,N,a\,\rangle$). As detailed in Ref. ​​[97], the relationship between $\hat{H}_{\mathrm{m-l}}$ and its form based on a single particle picture reveals a route to evaluate the coupling parameter V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b}. The evaluation of V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b} has been discussed in the Supporting Information of Ref. [​​68] and is briefly reviewed in Appendix A.

The molecule-lead coupling $\hat{H}_{\mathrm{m-l}}$ can also be written in terms of bath operators and system operators asGalperin and Nitzan (2005); Fainberg, Jouravlev, and Nitzan (2007); May and Kühn (2008a)

[TABLE]

where the definitions of system operators and bath operators are respectively given by

[TABLE]

and

[TABLE]

with the help of the decomposition of the coupling parameterPeskin (2017)

[TABLE]

$T^{\alpha}_{NaN-1b}$ refers to the probability amplitude of a molecular transition due to electrode $\alpha$, and M_{\alpha\mbox{\displaystyle\textbf{k}}\sigma} determines the spectral density of electrode $\alpha$.

II.2 Polaron Transformation

In order to simplify the subsequent derivation of master equations, we employ a polaron transformationNitzan (2006); Mahan (2008); May and Oliver (2011). We choose a representation of density matrix operator that is defined with respect to the electron-vibrational states $\displaystyle|\,N,a,\mu\,\rangle$. This electron-vibrational state is expressed as a tensor product of a many-electron state and an associated displaced harmonic oscillator state, i.e.,

[TABLE]

where the vibrational quantum number $\nu=0,1,2,3\dots$ and the corresponding energy level $E_{Na\nu}=E_{N,a}+(\nu+1/2)\hbar\omega_{\mathrm{vib}}$. Note that the vibrational operators $\hat{b}$ and $\hat{b}^{\dagger}$ in Eq. (II.1) are defined for \{\mbox{\displaystyle|,\nu\mbox{$\displaystyle\left(0\right)$},\rangle}\}, which is related to the states of a displaced harmonic oscillator \{\mbox{\displaystyle|,\nu\mbox{$\displaystyle\left(\lambda\right)$},\rangle}\} through \mbox{\displaystyle|,\nu\mbox{$\displaystyle\left(\lambda\right)$},\rangle}=\mbox{\displaystyle\mathbf{e}^{\lambda\mbox{$\displaystyle\left(\hat{b}^{\dagger}-\hat{b}\right)$}}}\mbox{\displaystyle|,\nu\mbox{$\displaystyle\left(0\right)$},\rangle}.Cahill and Glauber (1969)

After a polaron transformationKoch and von Oppen (2005); Mahan (2008); Hsu, Tsai, and Jin (2010a), the new molecular Hamiltonian $\hat{H}^{\prime}_{\mathrm{mol}}=\hat{U}\hat{H}_{\mathrm{mol}}\hat{U}^{-1}$ is given in terms of electron-vibrational states $\displaystyle|\,N,a,\mu\,\rangle$ as

[TABLE]

where the unitary operator $\hat{U}$ is

[TABLE]

The coupling of the molecule to the thermal bath in Eq. (17) is transformed similarly as

[TABLE]

with $\hat{B}^{\prime}_{\mathrm{m-th}}=\hat{B}_{\mathrm{m-th}}$ and

[TABLE]

Furthermore, $\hat{U}$ allows us to transform the coupling Hamiltonians $\hat{H}_{\mathrm{field}}$, $\hat{H}_{\mathrm{m-l}}$ and $\hat{H}_{\mathrm{SOC}}$ into $\hat{H}^{\prime}_{\mathrm{field}}$, $\hat{H}^{\prime}_{\mathrm{m-l}}$ and $\hat{H}^{\prime}_{\mathrm{SOC}}$ in the polaron frame by using the following relationship,

[TABLE]

The transformed molecule-lead coupling Hamiltonian then reads

[TABLE]

with $B^{\prime\pm}_{\alpha}=B^{\pm}_{\alpha}$ and

[TABLE]

where M_{\nu_{a}\nu_{b}}\mbox{\displaystyle\left(\lambda_{N-1,b}-\lambda_{N,a}\right)} absorbs the exponential factor introduced in Eq. (II.2) and refers to the vibrational overlap \displaystyle\langle\,\nu_{a}\mbox{\displaystyle\left(\lambda_{N,a}\right)}\,|\,\nu_{b}\mbox{\displaystyle\left(\lambda_{N-1,b}\right)}\,\rangle. M_{\nu_{a}\nu_{b}}\mbox{\displaystyle\left(\lambda\right)} for real $\lambda$ is then given byCahill and Glauber (1969); Koch et al. (2004)

[TABLE]

where \mathrm{sgn}\mbox{\displaystyle\left(x\right)} returns the sign of $x$, $\nu_{\mathrm{min}}=\mathrm{min}\{\nu_{a},\nu_{b}\}$, $\nu_{\mathrm{max}}=\mathrm{max}\{\nu_{a},\nu_{b}\}$ and \mathrm{L}^{n}_{m}\mbox{\displaystyle\left(x\right)} denotes the generalized Laguerre polynomial.

Similarly, one can derive the radiative coupling Hamiltonian $\hat{H}_{\mathrm{field}}$ and the SOC Hamiltonian $\hat{H}_{\mathrm{SOC}}$ in the polaron frame as

[TABLE]

and

[TABLE]

III Master Equation Approach

In this section, we will show that the electron transport characteristics of an irradiated molecular junction can be described using Pauli master equations (rate equations). The Pauli master equations built upon electron-vibrational states have been extensively discussedKoch et al. (2004); Koch and von Oppen (2005); Hsu, Tsai, and Jin (2010b). In most studies, the rate equations of incoherent sequential tunneling are phenomenologically constructed based on the Fermi golden rule in the framework of a T-matrix approach.Averin and Nazarov (1990, 1992); Esteve (1992); Schoeller and Schön (1994); Averin (1994); Matveev and Andreev (2002); Koch et al. (2004); Koch and von Oppen (2005); Walldorf, Jauho, and Kaasbjerg (2017) In our work, we derive the rate equations for irradiated incoherent sequential tunneling explicitly from the reduced density matrix approachBlum (1981); Cohen-Tannoudji, Dupont-Roc, and Grynberg (1992); Nitzan (2006); May and Oliver (2011) via the Redfield theoryRedfield (1965), incorporating the optical transitions, vibrational relaxations and singlet-triplet transitions.

First of all, the reduced density operator of the molecular system is defined as \hat{\rho}\mbox{\displaystyle\left(t\right)}=\mathrm{Tr}_{\mathrm{lead+th}}\{\hat{\rho}_{\mathrm{tot}}\}, where $\hat{\rho}_{\mathrm{tot}}$ is the density operator of the total system and the trace $\mathrm{Tr}_{\mathrm{lead+th}}\{\cdots\}$ averages over the states of the lead electrons and the thermal bath modes. The population of a state $\displaystyle|\,N,a,\nu\,\rangle$ is then defined as P_{N,a,\nu}\mbox{\displaystyle\left(t\right)}=\mbox{\displaystyle\langle,N,a,\nu,|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,N,a,\nu,\rangle}

We partition the total Hamiltonian $\hat{H}_{\mathrm{tot}}$ as in Eq. (1) and obtain the Liouville equation of the density operator $\hat{\rho}_{\mathrm{tot}}$ in the interaction picture as,

[TABLE]

where an operator in the interaction picture \hat{A}^{I}\mbox{\displaystyle\left(t\right)} relates to its Schrodinger picture counterpart \hat{A}\mbox{\displaystyle\left(t\right)} through \hat{A}^{I}\mbox{\displaystyle\left(t\right)}=\mbox{\displaystyle\mathbf{e}^{\frac{i}{\hbar}\hat{H}{0}t}}\hat{A}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle\mathbf{e}^{-\frac{i}{\hbar}\hat{H}{0}t}}.

Due to different coupling terms, the evolution of the density operator can be separated into

[TABLE]

where the four terms on the right hand side are respectively denoted as \frac{d\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}}{dt}|_{\mathrm{m-l}}, \frac{d\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}}{dt}|_{\mathrm{m-th}}, \frac{d\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}}{dt}|_{\mathrm{field}} and \frac{d\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}}{dt}|_{\mathrm{SOC}} in the following.

Based on Eq. (III), we will derive the rate equations of $\{P_{Na\nu}\}$ separately according to $\hat{H}^{I}_{\mathrm{m-l}}$, $\hat{H}^{I}_{\mathrm{m-th}}$, $\hat{H}^{I}_{\mathrm{field}}$ and $\hat{H}^{I}_{\mathrm{SOC}}$. Note that the rate equations due to $\hat{H}_{\mathrm{SOC}}$ are not derived explicitly. Instead, we treat the rates due to $\hat{H}_{\mathrm{SOC}}$ phenomenologically, because the relevant processes are insignificant compared to the other electronic transitions.

The molecule-lead rate equations $\frac{dP_{Na\nu}}{dt}|_{\mathrm{m-l}}$ and the vibrational relaxation rate equation $\frac{dP_{Na\nu}}{dt}|_{\mathrm{m-th}}$ will be derived through a 2nd order expansion of the Liouville equation using Redfield theory, whereas the derivation of the radiative rate equations $\frac{dP_{Na\nu}}{dt}|_{\mathrm{field}}$ does not require a perturbation expansion since the radiative coupling is semiclassically treated in $\hat{H}_{\mathrm{field}}$.

The rest of this section is organized as follows. In section III.1, we briefly review Redfield theory. In section III.2, III.3 and III.4, we respectively discuss the derivations of charge transfer transitions, vibrational relaxations and optical transitions. Section III.5 finalizes the derivation of the Pauli master equation via phenomenologically introducing the other relevant processes. Finally, section III.6 presents the current formula provided by the master equation approach.

III.1 Redfield Theory

We start from a component of the quantum master equation in Eq. (III), i.e.,

[TABLE]

where "coup" refers to the contributions from $\hat{H}_{\mathrm{m-l}}$ or $\hat{H}_{\mathrm{m-th}}$ by "m-l" or "m-th".

We then substitute the solution of Eq. (36), i.e., \hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}=\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t_{0}\right)}-\frac{i}{\hbar}\int_{t_{0}}^{t}du[\hat{V}^{I}\mbox{\displaystyle\left(u\right)}\,,\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(u\right)}], into Eq. (38), apply $\mathrm{Tr}_{\mathrm{lead+th}}\{\cdots\}$ on both sides, and employ the Born-Markov approximation. Finally we arrive at the following Redfield master equation,

[TABLE]

When the molecule is weakly coupled to the electrodes and the thermal bath, the electrons in the electrodes and the thermal bath modes are supposed to be relaxed significantly faster than all the other relevant processes, and the two baths are therefore assumed to be in their thermal equilibrium characterized by $\hat{\rho}_{\mathrm{lead}}$ and $\hat{\rho}_{\mathrm{th}}$ at all times. As a result, the Born-Markov approximation required by Redfield theory is valid in our problem. The Born approximation \hat{\rho}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}\approx\hat{\rho}\mbox{\displaystyle\left(t\right)}\hat{\rho}_{\mathrm{lead}}\hat{\rho}_{\mathrm{th}} is first invoked to decouple the dynamics of the molecular system from that of the lead electrons and the thermal bath. The non-local memory effect is then disregarded by replacing \hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(u\right)} with \hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}, which is the first Markov approximation. The second Markov approximation follows by increasing the upper limit of the integral to infinity after a variable transformation from $u$ to $u=t-u$.

In Eq. (III.1), it should be noted that the first order term \mathrm{Tr}_{\mathrm{coup}}\{\mbox{\displaystyle\left[\hat{H}{\mathrm{m-l}}\mbox{$\displaystyle\left(t\right)$},,\hat{\rho}^{I}{\mathrm{tot}}\mbox{$\displaystyle\left(t_{0}\right)$}\right]}\} has been eliminated by invoking an uncoupled initial state \hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t_{0}\right)}=\hat{\rho}^{I}\mbox{\displaystyle\left(t_{0}\right)}\hat{\rho}_{\mathrm{lead}}\hat{\rho}_{\mathrm{th}}. Besides, the 2nd order terms including the mix of $\hat{H}_{\mathrm{coup}}$ and other coupling Hamiltonians are excluded because they do not survive after the trace operation $\mathrm{Tr}_{\mathrm{lead+th}}\{\cdots\}$.

III.2 Charge Transfer Transition

In order to derive the charge transfer transitions, we first write down the transformed molecule-lead coupling Hamiltonian in the interaction picture as

[TABLE]

with

[TABLE]

and

[TABLE]

A substitution of the $\hat{H}^{\prime I}_{\mathrm{m-l}}$ in Eq. (40) into Eq. (III.1) results in the following master equation,

[TABLE]

where two types of correlation functions for electrode electrons are defined as

[TABLE]

and

[TABLE]

The spectra density J_{\alpha}\mbox{\displaystyle\left(\omega\right)} of the electrons in electrode $\alpha$ is specified by

[TABLE]

and the Fermi distribution function f_{\beta}\mbox{\displaystyle\left(\epsilon,\mu\right)} is given by

[TABLE]

where $\mu$ is the electron chemical potential.

In order to derive the rate equations, we evaluate the diagonal matrix element with respect to electron-vibrational states \{\mbox{\displaystyle|,N,a,\nu,\rangle}\} on both sides of Eq. (III.2) and then arrive at

[TABLE]

where $\Re$ refers to the real part.

We next take the evaluation of the first term in Eq. (III.2) as an example, i.e.,

[TABLE]

which indicates that the evolution of \mbox{\displaystyle\langle,N,a,\nu,|}\hat{\rho}^{I}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,N,a,\nu,\rangle} (the population $P_{Na\nu}$) depends on the off-diagonal density matrix elements (coherence)Scholes et al. (2018).

We assume that the dynamics of the molecular system is resolved over a time step $\Delta t$ during which the reduced density matrix \hat{\rho}^{I}\mbox{\displaystyle\left(t\right)} does not vary significantly, then the terms that satisfy \Delta t\gg\mbox{\displaystyle\left|E_{N+1,a_{2}}-E_{N+1,a_{1}}+\mbox{$\displaystyle\left(\nu_{2}-\nu_{1}\right)$}\hbar\omega\right|}^{-1} does not contribute to the propagation of \hat{\rho}^{I}\mbox{\displaystyle\left(t\right)} due to fast oscillations. As a result, only secular terms satisfying $E_{N+1,a_{2}}-E_{N+1,a_{1}}+(\nu_{2}-\nu)\hbar\omega=0$ are kept. The consequence of the secular approximation differs for molecular systems with and without degenerate states. When there is no degeneracy in the system, the corresponding exponential factor with local time dependence can be replaced by Kronecker functions, e.g.,

[TABLE]

which directly decouples the evolution of the population terms from the coherence terms of \hat{\rho}^{I}\mbox{\displaystyle\left(t\right)}. However, when the molecular system contains degenerate states, such as the system considered in this paper, terms satisfying $\nu_{1}=\nu_{2}$ and $E_{N+1,a_{2}}=E_{N+1,a_{1}}$ are still left after the secular approximation has been applied, wherein the coherence between degenerate electronic states is involved in the evolution of the populations. However, considering that the electronic dephasing is the fastest process in the system, the coherence terms can be neglected when studying steady-state transport characteristics. Consequently, we can drop all the terms containing the off-diagonal density matrix elements in Eq. (III.2) regardless of the presence of degeneracy. As a result, the evaluation of the matrix elements gives us

[TABLE]

where the diagonal matrix element has been replaced with the population term due to the relationship

[TABLE]

We can use the same procedure as above to evaluate the other terms in Eq. (III.2). Finally, we obtain the following molecule-lead rate equations

[TABLE]

The rates of electron transport and hole transport processes are thus given by

[TABLE]

in which the prefactors read

[TABLE]

where $\bar{J}_{\alpha}$ denotes a constant spectral density of electrode $\alpha$ in the wide band limit. After evaluating $T^{\alpha}_{N\pm 1b,N,a}$ as detailed in Appendix A, we can rewrite the prefactor $\gamma^{\alpha}_{N-1b\nu^{\prime},Na\nu}$ as

[TABLE]

where \Gamma_{\alpha}=\frac{2\pi}{\hbar^{2}}\mbox{\displaystyle\left|\zeta_{\alpha}\right|}\bar{J}_{\alpha} is the characteristic rate of charge transfer transition associated with electrode $\alpha$, and the dimensionless coupling $\nu_{N\pm 1b,Na}$ of the transition \mbox{\displaystyle|,N\pm 1,b,\rangle}\leftrightarrow\mbox{\displaystyle|,N,a,\rangle} has been given in Table S1 in Ref. [​​68] for all relevant transitions considered in this paper. The meaning of $\zeta_{\alpha}$ and the evaluation of $\nu_{N\pm 1b,Na}$ are also discussed in the Appendix A.

III.3 Vibrational Relaxation

The derivation of the vibrational relaxation rate equations follows the same procedure as the above section. In the interaction picture, the transformed coupling Hamiltonian $\hat{H}^{\prime}_{\mathrm{m-th}}$ reads

[TABLE]

where

[TABLE]

and

[TABLE]

For simplicity, we abbreviate $\displaystyle|\,N,a\,\rangle$ as $\displaystyle|\,A\,\rangle$ here. We then plug Eq. (58) into Eq. (III.1) and obtain the following master equation,

[TABLE]

where the thermal bath correlation function is defined as

[TABLE]

The spectral density J_{\mathrm{th}}\mbox{\displaystyle\left(\omega\right)} of the thermal bath is defined as

[TABLE]

and the Bose-Einstein distribution function n_{\beta}\mbox{\displaystyle\left(\omega\right)} is given by

[TABLE]

A substitution of S^{\prime I}_{\mathrm{m-th}}\mbox{\displaystyle\left(t\right)} in Eq. (III.3) into Eq. (III.3) leads us to

[TABLE]

where rotating wave approximation has been applied, and $\hat{\rho}_{AA^{\prime}}=\langle A|\hat{\rho}|A^{\prime}\rangle$ refers to an electronic block of the reduced density matrix. Eq. (III.3) suggests that $\hat{H}_{\mathrm{m-th}}$ not only leads to the vibrational relaxations within the same electronic state manifold but also results in the dephasing of coherence between states associated with different electronic state manifold.

The vibrational relaxation rate equations $\frac{dP_{N,a,\nu}}{dt}|_{\mathrm{m-th}}$ is derived by taking the diagonal matrix element of Eq. (III.3), which gives

[TABLE]

where the rates of vibrational relaxation is

[TABLE]

Here, \gamma_{\mathrm{p}}=\frac{J_{\mathrm{th}}\mbox{\displaystyle\left(\omega_{\mathrm{vib}}\right)}}{m_{\mathrm{vib}}\omega_{\mathrm{vib}}} is the characteristics rate of vibrational relaxation. Note that the corresponding characteristic time scale $\tau_{\mathrm{p}}=\gamma_{\mathrm{p}}^{-1}$ refers to the life time of the first excited vibrational state.

III.4 Optical Transition

In the interaction picture, the transformed radiative coupling Hamiltonian reads

[TABLE]

Since the light-matter interaction is considered semi-classically, we derive the radiative rate equations without employing a $2$nd order expansion.

We start from \frac{d}{dt}\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}|_{\mathrm{field}}=-\frac{i}{\hbar}[\hat{H}^{I}_{\mathrm{field}},\hat{\rho}^{I}_{\mathrm{tot}}\mbox{\displaystyle\left(t\right)}], apply the Born approximation \hat{\rho}_{\mathrm{tot}}=\hat{\rho}\mbox{\displaystyle\left(t\right)}\hat{\rho}_{\mathrm{lead}}\hat{\rho}_{\mathrm{th}} and the trace operation $\mathrm{Tr}_{\mathrm{lead+th}}\{\cdots\}$. As a result, we have the following equation of motion of the reduced density operator \hat{\rho}^{I}\mbox{\displaystyle\left(t\right)},

[TABLE]

By substituting Eq. (III.4) into Eq. (69) and taking the diagonal matrix element with respect to $\displaystyle|\,0,S_{0},\nu_{0}\,\rangle$, one can derive the rate equation of $P_{0,S_{0},\nu_{0}}$,

[TABLE]

where $\mathrm{c.c.}$ stands for complex conjugate. In order to solve the coherence term \mbox{\displaystyle\langle,0,S_{1},\nu,|}\hat{\rho}^{I}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{0},\nu_{0},\rangle}, we turn to the following equation of motion in the Schrödinger picture,

[TABLE]

where $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ describes the pure dephasing of the electronic coherence which originates from the other coupling Hamiltonians, i.e., $\hat{H}_{\mathrm{m-l}}$, $\hat{H}_{\mathrm{m-th}}$. An explicit evaluation in Appendix B shows that $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ is a sum over the rates of all the processes that start from $\displaystyle|\,0,S_{0},\nu_{0}\,\rangle$ and $\displaystyle|\,0,S_{1},\nu_{1}\,\rangle$ excluding optical transitions, e.g., charge transfer transitions and vibrational relaxations.

When only the steady-state dynamics is concerned, \mbox{\displaystyle\langle,0,S_{1},\nu^{\prime},|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{1},\nu_{1},\rangle} and \mbox{\displaystyle\langle,0,S_{0},\nu_{0},|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{0},\nu,\rangle} in Eq. (III.4) are time-independent. In this case, according to the time-periodicity of E\mbox{\displaystyle\left(t\right)}, one can take a trial solution of \mbox{\displaystyle\langle,0,S_{0},\nu_{0},|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{1},\nu_{1},\rangle},

[TABLE]

\rho_{0S_{0}\nu_{0},0S_{1}\nu_{1}}\mbox{\displaystyle\left(\omega\right)} is then solved as

[TABLE]

where we have dropped the fast oscillating terms containing $\displaystyle\mathbf{e}^{\pm 2i\omega}$ (rotating wave approximation). After substituting the solution Eq. (72) and (III.4) into Eq. (III.4), we obtain the rate equation of $P_{0,S_{0},\nu_{0}}$ in the following,

[TABLE]

The coherence terms \mbox{\displaystyle\langle,0,S_{a},\nu,|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{a},\nu^{\prime},\rangle} in Eq. (III.4) can be dropped by transforming Eq. (III.4) back to the interaction picture and applying the secular approximation. We finally obtain the radiative rate equation of $\displaystyle|\,0,S_{0},\nu_{0}\,\rangle$, i.e.,

[TABLE]

with the rates of stimulated optical transitions $k_{0,S_{0},\nu_{1}\leftrightarrow 0,S_{1},\nu_{1}}^{\mathrm{field}}$ given by

[TABLE]

where $\Im$ refers to the imaginary part.

The rate equation of the other singlet state $\displaystyle|\,0,S_{1},\nu_{1}\,\rangle$ could be obtained similarly as

[TABLE]

III.5 Pauli Master Equation

The final form of the Pauli master equations is achieved by phenomenologically introducing the rates of singlet-triplet transitions and the rate of spontaneous emission. The singlet-triplet transitions include the processes of intersystem crossing \mbox{\displaystyle|,0,S_{1},\nu,\rangle}\rightarrow\mbox{\displaystyle|,0,T^{m}{1},\nu^{\prime},\rangle} and phosphorescence \mbox{\displaystyle|,0,T^{m}{1},\nu,\rangle}\rightarrow\mbox{\displaystyle|,0,S_{0},\nu^{\prime},\rangle}. We collect all the rates discussed above and arrive at the following form of Pauli master equations,

[TABLE]

where $k^{\mathrm{SOC}}_{N,a,\nu\leftarrow N^{\prime},a^{\prime},\nu^{\prime}}$ refers to the rate due to SOC and $k^{\mathrm{spon}}_{0,S_{0},\nu\leftarrow 0,S_{1},\nu^{\prime}}$ refers to the rate of spontaneous emission.

III.6 Current Formula

The electric current through the molecular junction is defined by the flow of electron through electrode $\alpha$, i.e.,

[TABLE]

where the trace $\mathrm{Tr}\{\cdots\}$ includes the average over the electrons in the electrodes, the thermal bath modes and the electron-vibrational states of the molecular system. Following the same procedure for deriving the molecule-lead rate equations, we obtain the following current formula,

[TABLE]

Since we study the steady-state current, the population $P_{N,a,v}$ takes the stationary solution of the Pauli master equation.

The current formula Eq.​​​​​ (80) suggests that the current through the molecular junction is determined by both the rates of the charge transfer transitions and the steady-state population of the involved states, wherein the population is determined by the rates of relevant electronic transitions. It is noted that the source-drain bias $V_{\mathrm{SD}}$ and the gate voltage $V_{\mathrm{G}}$ control the rates of charge transfer transitions in Eqs. (54) and (55) by virtue of the determination of the Fermi distribution function. In the low temperature limit, the Fermi distribution function becomes a step function. In this case, the rates of transitions between a neutral state $\displaystyle|\,0,a,\nu_{a}\,\rangle$ and a charged state $\displaystyle|\,N,b,\nu_{b}\,\rangle$ ($N=\pm 1$) can be clearly expressed as a function of $V_{\mathrm{SD}}$ and $V_{\mathrm{G}}$,

[TABLE]

and

[TABLE]

The step function \theta\mbox{\displaystyle\left(x\right)} takes a value of $1$ for $x>0$ and a value of [math] otherwise. The renormalized state energy associated with electrode $\alpha$ is defined by

[TABLE]

where $\zeta_{\mathrm{S}}=1$ and $\zeta_{\mathrm{D}}=-1$. Eq. (83) combines the energy level of $\displaystyle|\,N,a,\nu\,\rangle$ with the Fermi energy on electrode $\alpha$, i.e., $\mu_{0}+\zeta_{\alpha}eV_{\mathrm{SD}}/2$, and the energy shift caused by gate voltage, i.e., $-NeV_{\mathrm{G}}$. Eqs. (81) and (82) clearly reveal that the charge transfer transitions \mbox{\displaystyle|,0,a,\nu_{a},\rangle}\leftrightarrow\mbox{\displaystyle|,N,b,\nu_{b},\rangle} ($N=\pm 1$) are explicitly determined by the energy level alignments between $E_{0,a,\nu}$ and \tilde{E}^{\alpha}_{N,a,\nu}\mbox{\displaystyle\left(V_{\mathrm{G}},V_{\mathrm{SD}}\right)}, which provides us with a handy tool for identifying the dominant transport channels.

IV Parameters

We compute the electron transport characteristics in the succeeding section using the system parameters compiled in Table 1. As shown in Figure 2, we consider a model system with a specific configuration of PESs. In this configuration, the electronic transitions between neutral states are not coupled to the reaction coordinate, whereas the PESs of cation states and anion states are displaced oppositely with respect to the PESs of neutral states by the same amount. The displacement between the PESs is characterized by the e-p coupling parameter $\lambda$. The energies of the involved many-electron states \mbox{\displaystyle|,N,a,\rangle}\in{\cal F}, which corresponds to the local minimum of the associated PESs as shown in Figure 2, are determined by the ab initio calculations of zinc phthalocyanine at the level of TDDFT/B3LYP/6-311g(d,p) using Gaussian 16.Frisch et al. (2016) The charged state energies are corrected for the image charge effect, see Supporting Information of Ref.  [​​68] for details. A value of $\mu_{0}=-5.3$ eV is adopted for the chemical potential of Au(111) facet throughout all the computed transport characteristics. A crossover from weak to strong e-p couplings will be explored using $\lambda=0.25,0.5,1.0,2.0$. In order to comply with the weak molecule-lead coupling limit, we adopt a value of $0.0001$ eV for $\hbar\Gamma$ with symmetric molecule-lead couplings $\Gamma_{\mathrm{S}}=\Gamma_{\mathrm{D}}=\Gamma$. The energy of vibrational quanta is chosen as $0.2$ eV and $0.02$ eV. Moreover, we investigate electron transport characteristics in the low temperature regime ($\beta=0.05\hbar\omega_{\mathrm{vib}}$) in order to resolve the Frank-Condon allowed charge transfer transitions in the stability diagram and the current-voltage characteristics.

Since there is no displacement between electronic manifold $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,0,S_{1}\,\rangle$, a value of $10^{12}$ $\mathrm{s}^{-1}$ is chosen for $k^{\mathrm{field}}_{0,S_{0},\nu\leftrightarrow 0,S_{1},\nu}$ such that $k^{\mathrm{field}}\Gamma^{-1}\gtrsim 1$, where the rate of optical transition is hereinafter referred to as $k^{\mathrm{field}}$. Assuming the optical excitation is in resonance with the electronic transition \mbox{\displaystyle|,0,S_{0},\nu,\rangle}\leftrightarrow\mbox{\displaystyle|,0,S_{1},\nu,\rangle}, the respective electric field strength takes the value of $100$ $\mathrm{kV}/\mathrm{cm}$ if the transition dipole moment is $5$ Debye. A spontaneous emission rate of $10^{8}$ $\mathrm{s}^{-1}$ estimated through 4\alpha n\omega_{10}^{3}|\mbox{\displaystyle\textbf{d}}_{01}|^{2}/(3c^{2}) is also introduced phenomenologically, where $\alpha$ is the fine structure constant, $\omega_{10}=(E_{0,S_{1}}-E_{0,S_{0}})/\hbar$, $n$ is the refractive index and $c$ is the speed of light in vacuum. The rates of singlet-triplet transitions are then chosen for typical values, i.e., $10^{6}$ $\mathrm{s}^{-1}$ for intersystem crossing processes \mbox{\displaystyle|,0,S_{1},\rangle}\rightarrow\mbox{\displaystyle|,0,T^{m}{1},\rangle} and $10^{3}$ $\mathrm{s}^{-1}$ for phosphorescent processes \mbox{\displaystyle|,0,T^{m}{1},\rangle}\rightarrow\mbox{\displaystyle|,0,S_{0},\rangle}.

V Results and Discussion

In this section, we investigate the influence of e-p coupling and vibrational relaxation on the transport characteristics under an optical field. Light-driven electron transport is explored from weak to strong e-p coupling regimes. For simplicity, the effect of vibrational relaxation is discussed in two limits: the equilibrated phonon and the unequilibrated phonon. We analyze the transport characteristics from the following aspects. In section V.2, we study the charge stability diagram (the conductance spectra) of the irradiated molecular junction in a broad range of source-drain bias voltage and gate voltage. Charge stability diagrams have been extensively studied in nanoscale electron transport and they can provide rich information about electronic structures. Moreover, considering that low-bias current-voltage characteristics is accessible in most experiments, we focus on low-bias transport characteristics in section V.3 and V.4. In section V.3, we control the gate voltage, compute the low-bias current at $V_{\mathrm{SD}}=0.1$ V, and compare the situations of high- and low-frequency vibrational modes in a range of gate voltage $V_{\mathrm{G}}$ that covers all relevant energy level alignment schemes between the charged states $\displaystyle|\,-1,D^{\alpha}_{0}\,\rangle$ and the neutral states. In section V.4, focusing on four representative schemes of energy level alignments, we compute the current-voltage characteristics and derive the analytical current formula. This offers a quantitative point of view on the light-driven transport when the charge transfer transitions are coupled to intramolecular vibrations.

V.1 Equilibrated and Unequilibrated Phonon regimes

These two extreme limits of vibrational relaxation, i.e., $\gamma_{\mathrm{p}}\gg\Gamma$ (equilibrated phonon regime) and $\gamma_{\mathrm{p}}\ll\Gamma$ (unequilibrated phonon regime), offer us a clear picture to understand the role of the vibrational relaxation on the light-driven transport characteristics.

In the equilibrated phonon regime, the vibrational relaxation is much faster than all the other transitions about the molecule. In other words, the vibrational state distribution of each electronic state manifold instantaneously relax into its thermal equilibrium upon any electronic transition. As a result, the vibrational relaxation part of the Pauli master equation Eq. (III.5), i.e., $\frac{dP_{N,a,\nu}}{dt}|_{\mathrm{m-th}}$, reduces into

[TABLE]

where P^{\mathrm{eq}}_{\nu}=\mbox{\displaystyle\mathbf{e}^{\nu\beta\hbar\omega}}/\sum_{\nu^{\prime}=0}^{\infty}\mbox{\displaystyle\mathbf{e}^{\nu^{\prime}\beta\hbar\omega}} describes the equilibrium vibrational distribution. The role of Eq. (84) on the Pauli master equation Eq. (III.5) is to force the vibrational state distribution to be kept at the thermal equilibrium distribution at all times. This matter of fact enables us to derive analytial solutions for transport characteristics in section V.4.

In the unequilibrated phonon regime, on the contrary, we assume that the vibrational relaxation is slower than any other processes about the molecule, which simply modifies the Pauli master equation Eq. (III.5) by dropping the terms belonging to $\frac{dP_{Na\nu}}{dt}|_{\mathrm{m-th}}$. In other words, the non-equilibrium vibrational excitations brought by other electronic transitions are completely preserved during the time scale of the relevant processes. The role of the non-equilibrium vibrational excitations are then investigated in this regime.

V.2 Charge Stability Diagram: Effect of e-p Coupling on Photoinduced Coulomb Diamond

We explore the influence of e-p coupling ($\lambda$) on the charge stability diagram, i.e., a plot of differential conductance $dI/dV$ versus $V_{\mathrm{SD}}$ and $V_{\mathrm{G}}$, in the equilibrated phonon regime and the unequilibrated phonon regime Koch and von Oppen (2005). We discuss the equilibrated phonon regime in details. For the unequilibrated phonon regime, we only discuss the difference from the equilibrated phonon regime.

We first investigate the equilibrated phonon regime, which is characterized by the equilibrium vibrational distribution. As presented in Figure 3, we compute the conductance spectra for $\hbar\omega_{\mathrm{vib}}=0.2$ eV and $\lambda=0.25,0.5,1.0,2.0$ under a field-off condition $k^{\mathrm{field}}=0$ and a field-on condition $k^{\mathrm{field}}=10^{12}$ $\mathrm{s}^{-1}$.

When the radiation is off, the conventional Coulomb blockade diamonds A, B and C, respectively corresponding to $\displaystyle|\,-1,D^{\pm 1/2}_{0}\,\rangle$, $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,1,D^{\pm 1/2}_{0}\,\rangle$, remain intact for e-p couplings from $\lambda=0.25$ to $\lambda=1.0$, as shown in Figure 3a-c. The diamond B, corresponding to $\displaystyle|\,0,S_{0}\,\rangle$, exhibits the width $E_{-1,D^{\sigma}_{0}}+E_{1,D^{\sigma}_{0}}-2E_{0,S_{0}}$ and the height-to-width ratio $2:1$. From Figure 3a to Figure 3c, as the e-p coupling $\lambda$ increases, conductance lines spaced by $\hbar\omega_{\mathrm{vib}}$ along $V_{\mathrm{G}}$ axis and $2\hbar\omega_{\mathrm{vib}}$ along $V_{\mathrm{SD}}$ axis become pronounced. The equally spaced conductance lines arise from the step-wise activation of charge transport channels. In Figure 3d, when $\lambda$ becomes $2.0$, we observe a clear Frack-Condon blockade pattern, which has been reported both theoreticallyKoch and von Oppen (2005) and experimentallyBurzuri et al. (2014). The two intersections in between the three Coulomb diamonds break because the low-bias conduction is suppressed by the Frank-Condon factor (M_{\nu_{1}\nu_{2}}\mbox{\displaystyle\left(\lambda\right)} in Eqs. (54) and (55)), i.e., the transitions between low-lying vibrational states decrease exponentially in the regime of strong e-p coupling ($\lambda>1$).

When the radiation is on, the anomalous Coulomb blockade pattern is clearly identified in the weak e-p coupling regime ($\lambda<1.0$) from Figure 3e to Figure 3g. In the zero e-p coupling limit, the widths of the diamonds D, E and F are $E_{0,T^{m}_{1}}-E_{0,S_{0}}$, $E_{-1,D^{\sigma}_{0}}+E_{1,D^{\sigma}_{0}}-2E_{0,T^{m}_{1}}$ and $E_{0,T^{m}_{1}}-E_{0,S_{0}}$, respectively. As the e-p coupling increases, the diamonds D and F (corresponding to partial charged states) shrink by steps of \Delta V_{\mathrm{G}}=\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|} ,while the diamond E (corresponding to the triply degenerate states $\displaystyle|\,0,T^{m}_{1}\,\rangle$) remains invariant. The conductance lines in the diamonds D and F result from the additional transport channels activated by e-p coupling, while the robustness of the diamond E with respect to the e-p coupling is attributed to the origin of the diamond E, i.e., the presence of the triplet states.

In the limit of unequilibrated phonon regime, the field-off and field-on charge stability diagrams (conductance spectra) are computed for $\lambda=0.25,0.5,1.0,2.0$, see Figure 4. A comparison between Figure 3 and Figure 4 shows that the non-equilibrium vibrational exicitations result in more reduction in the sizes of diamonds D and F at the same e-p coupling. At a low bias, the equilibrated and unequilibrated phonon regimes exhibit the same trend in response to e-p coupling and optical excitation.

V.3 Current vs Gate Voltage: Role of Vibrational Frequency and Vibrational Relaxation

In most cases, the transport characteristics of a molecular junction is only available at a low bias voltage due to the instability caused by the high-bias electric field. Therefore, in this section, we concentrate on the field-on current at a low bias, i.e., $V_{\mathrm{SD}}=0.1$ V, and explore the current response via changing the gate voltage $V_{\mathrm{G}}$. Similar to the previous section, we carry out the calculation for both the equilibrated phonon regime and the unequilibrated phonon regime. Furthermore, in additional to the situation with high frequency mode $\hbar\omega_{\mathrm{vib}}=0.2$ eV, we also consider a situation in which a low frequency vibrational mode ($\hbar\omega_{\mathrm{vib}}=0.02$ eV) dominates. The results are presented in Figure 5. In view of the symmetry in the charge stability diagrams, we only present the current within a range of gate voltages that covers all possible energy level alignment schemes between \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)} and neutral states. For simplicity, we hereafter denote the gate voltage corresponding to \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0} as $V^{\mathrm{S0}}_{\mathrm{G}}$ and the gate voltage corresponding to \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,T^{m}_{1},0} as $V^{\mathrm{T1}}_{\mathrm{G}}$.

Generally, from Figure 5a to 5d, we identify the plateaus of current centered at the gate voltages that correspond to the energy level alignments \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0}+n\hbar\omega_{\mathrm{vib}} ($n$ is non-negative integer) and \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,T^{m}_{1},0}. At zero e-p coupling, the electronic transitions \mbox{\displaystyle|,-1,D^{\sigma}{0},\nu,\rangle}\leftrightarrow\mbox{\displaystyle|,0,S{0},\nu^{\prime},\rangle} and \mbox{\displaystyle|,-1,D^{\sigma}{0},\nu,\rangle}\leftrightarrow\mbox{\displaystyle|,0,T^{m}{1},\nu^{\prime},\rangle} are allowed for $\nu=\nu^{\prime}$, so we observe the current plateaus only at $V^{\mathrm{S0}}_{\mathrm{G}}$ and $V^{\mathrm{T1}}_{\mathrm{G}}$. As the e-p coupling increases, the current at these two plateaus decrease monotonically and additional photoinduced current plateaus show up between $V^{\mathrm{S0}}_{\mathrm{G}}$ and $V^{\mathrm{T1}}_{\mathrm{G}}$, because the e-p coupling suppresses the diagonal vibrational transitions, i.e., the electronic transitions with $\nu=\nu^{\prime}$, and invokes the off-diagonal vibrational transitions, i.e., the electronic transitions with $\nu\neq\nu^{\prime}$.

Considering a high-frequency mode, e.g., $\hbar\omega_{\mathrm{vib}}=0.2$ eV, we observe four isolated current plateaus between $V^{\mathrm{S0}}_{\mathrm{G}}$ and $V^{\mathrm{T1}}_{\mathrm{G}}$ in Figure 5a and 5b. The current plateaus are isolated from each other because \hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}>V_{\mathrm{SD}}. A comparison between Figure 5a and 5b demonstrates that, at the same e-p coupling, the non-equilibrium vibrational population in the unequilibrated phonon regime leads to more pronounced current plateaus within $[V^{\mathrm{S0}}_{\mathrm{G}},V^{\mathrm{T1}}_{\mathrm{G}}]$.

In order to find out the origin of the e-p coupling induced current plateaus within $[V^{\mathrm{S0}}_{\mathrm{G}},V^{\mathrm{T1}}_{\mathrm{G}}]$, we consider a low-frequency mode, i.e., $\hbar\omega_{\mathrm{vib}}=0.02$ eV. In this case, the e-p coupling induced current plateaus are observed only in the vicinity of $V^{\mathrm{S0}}_{\mathrm{G}}$ and $V^{\mathrm{T1}}_{\mathrm{G}}$, because \hbar\omega_{\mathrm{vib}}\ll\mbox{\displaystyle\left|V^{\mathrm{T1}}{\mathrm{G}}-V^{\mathrm{S0}}{\mathrm{G}}\right|}. As a result, the role of the e-p coupling behaves like a broadening of the current plateaus at $V^{\mathrm{S0}}_{\mathrm{G}}$ and $V^{\mathrm{T1}}_{\mathrm{G}}$. In both extreme limits of vibrational relaxation, the additional Franck-Condon allowed transport channels is activated when $\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}$ is aligned with vibrational excitations of the singlet ground state. This broadens the current plateau at $V^{\mathrm{S0}}_{\mathrm{G}}$ toward higher gate voltages. In contrast, the e-p coupling induced current plateaus due to \mbox{\displaystyle|,-1,D^{\sigma}{0},\nu,\rangle}\leftrightarrow\mbox{\displaystyle|,0,T^{m}{1},\nu^{\prime},\rangle}, i.e., the broadening of the current plateau at $V^{\mathrm{T1}}_{\mathrm{G}}$, requires the nonequilibrium vibrational population and shows up only in the unequilibrated phonon regime.

V.4 Current vs Source-Drain Voltage: Energy Level Alignment and Analytic Analysis

In our previous workFu et al. (2018), we obtained the analytical solutions to the current-voltage characteristics without considering e-p coupling in the cases of several specific energy level alignments. These analytical solutions provide an in-depth and quantitative understanding on the transport characteristics of an irradiated molecular junction.

When the e-p coupling is significant, it is challenging to derive general analytical solutions, because the Pauli master equations are constructed within an enormous Fock space spanned by the electron-vibrational states. In order to analyze the effect of e-p coupling, for simplicity, we focus on the equilibrated phonon regime, which allows us to work on a smaller Fock space that is spanned by the electronic states. The equilibrated phonon regime enables the analytical solutions for two reasons. First, this regime excludes the vibrational relaxations, thus allowing us to consider only the transitions between different electronic manifolds, i.e., \mbox{\displaystyle|,N,a,\rangle}\rightarrow\mbox{\displaystyle|,N^{\prime},a^{\prime},\rangle}. Second, the equilibrium vibrational state distribution on each electronic state manifold is dominated by the vibrational ground state, because $\hbar\omega_{\mathrm{vib}}\gg k_{B}T$ is valid here.

Now we reconstruct the rate equations built upon the Fock space \{\mbox{\displaystyle|,N,a,\rangle}\} in the equilibrated phonon regime. Considering that $\hbar\omega_{\mathrm{vib}}\gg k_{B}T$, one can equate the population of an electronic state $\displaystyle|\,N,a\,\rangle$ to the population of its associated lowest lying vibrational state $\displaystyle|\,N,a,0\,\rangle$, i.e., $P_{N,a}=\sum_{\nu}P_{N,a,\nu}=P_{N,a,0}$. Therefore, the effective rate of an electronic transition \mbox{\displaystyle|,N,a,\rangle}\rightarrow\mbox{\displaystyle|,N^{\prime},a^{\prime},\rangle} is the sum of the rates of all outgoing processes from $\displaystyle|\,N,a,0\,\rangle$, i.e.,

[TABLE]

In this way, the effective rate of the charge transfer transitions can be expressed as

[TABLE]

where

[TABLE]

refers the ratio of the effective rate k^{\alpha\mbox{\displaystyle\left(\mathrm{eff}\right)}}_{N^{\prime},a^{\prime}\leftarrow N,a} to its zero e-p coupling limit. Note that $\Lambda^{\alpha}_{N^{\prime},a^{\prime}\leftarrow N,a}$ could be alternatively expressed as

[TABLE]

where

[TABLE]

and

[TABLE]

with the floor function $\lfloor x\rfloor$ returning the largest integer less than $x$. Eq. (89) indicates that $\Lambda^{\alpha}_{N^{\prime},a^{\prime}\leftarrow N,a}$ ranges from 0 to 1.0. Besides, $\Lambda^{\alpha}_{N^{\prime},a^{\prime}\leftarrow N,a}$ reaches its upper limit when \tilde{E}^{\alpha}_{Na0}\mbox{\displaystyle\left(V_{\mathrm{G}},V_{\mathrm{SD}}\right)}\gg\tilde{E}^{\alpha}_{N^{\prime}a^{\prime}0}\mbox{\displaystyle\left(V_{\mathrm{G}},V_{\mathrm{SD}}\right)} due to the relationship

[TABLE]

In order to better understand the effect of energy level alignment on light-driven electron transport when coupled to vibrations, we analyze the analytical solutions of the current-voltage characteristics and the corresponding numerical calculations in the four representative cases specified by the following energy alignments: (1) \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0}, (2) E_{0,S_{0},0}<\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}<E_{0,T^{m}_{1},0}, (3) \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,T^{m}_{1},0} and (4) \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}>E_{0,T^{m}_{1},0}.

V.4.1 \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0}

In this situation, the relevant electronic state manifolds are $\displaystyle|\,-1,D^{\uparrow}\,\rangle$, $\displaystyle|\,-1,D^{\downarrow}\,\rangle$, $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,0,S_{1}\,\rangle$. When the e-p coupling is zero, the energy level alignment \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0} leads to a field-off electric current with a zero onset voltage, similar to the resonant tunneling described by the Landauer theory. Compared to the field-off transport characteristics, the field-on electron transport in this case exhibits a suppression in the current, hereinafter referred to as photoresistance.

We now move on to the cases with e-p couplings. One can construct the Pauli master equations based on the electronic states using the effective rates given in Eq. (86). The rate equations in a matrix form is

[TABLE]

where the occupation probabilities in the vector form are

[TABLE]

Moreover, based on the energy level alignment in Figure 7a, the rate matrix W is constructed as

[TABLE]

For simplicity, we hereinafter denote $\Lambda^{\alpha}_{0,S_{0}\leftrightarrow-1,D^{\uparrow/\downarrow}_{0}}$ as $\Lambda^{\alpha}_{0}$, $\Lambda^{\alpha}_{0,S_{1}\leftrightarrow-1,D^{\uparrow/\downarrow}_{0}}$ as $\Lambda^{\alpha}_{1}$ and neglect the insignificant processes, i.e., the singlet-triplet transitions and the spontaneous emission.

By solving the steady-state dynamics from \mbox{\displaystyle\textbf{W}}\cdot\mbox{\displaystyle\textbf{P}}=0, we obtain an analytical expression of the stationary current as

[TABLE]

where we employ the abbreviations $\Lambda_{0/1}=\Lambda^{\mathrm{S}}_{0/1}+\Lambda^{\mathrm{D}}_{0/1}$. Since $I_{\mathrm{S}}+I_{\mathrm{D}}=0$ is valid in stationary transport, we only refer to $I_{\mathrm{S}}$ below. It is worth noting that, when $\lambda\rightarrow 0$, the analytical expression Eq. (99) is reduced to \mbox{\displaystyle\left|e\right|}\Gamma\mbox{\displaystyle\left(3k^{\mathrm{field}}+4\Gamma\right)}/\mbox{\displaystyle\left(6k^{\mathrm{field}}+6\Gamma\right)}, which agrees with our previous workFu et al. (2018).

Since $\Lambda^{\mathrm{S}}_{0}=\Lambda^{\mathrm{D}}_{0}=0.5\Lambda_{0}$ in the situation \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0}, the analytical expression of electric current, Eq. (99), could be further reduced to

[TABLE]

which clearly reveals that the field-off current (Recall $k^{\mathrm{field}}=0$) is

[TABLE]

and the amount of current suppression led by optical excitation is

[TABLE]

In Figure 6a, the photoresistivity is clearly identified in all explored e-p coupled regimes. Obviously, the current decreases as the e-p coupling increases, which is also clearly described in Eq. (99).

Furthermore, the stepwise dependence of $\Lambda^{\alpha}_{0/1}$ on $V_{\mathrm{SD}}$ and $V_{\mathrm{G}}$ results in the steplike current-voltage characteristics with width of \Delta V_{\mathrm{SD}}=2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|} and the equally spaced conductance lines in the charge stability diagrams. In other words, from the aspects of the energy level alignment, the current changes whenever \tilde{E}^{\alpha}_{\pm 1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},V_{\mathrm{SD}}\right)} goes across an electron-vibrational state $\displaystyle|\,0,S_{0},\nu\,\rangle$. Therefore, a change of $\hbar\omega_{\mathrm{vib}}$ in the energy level alignment corresponds to either a \Delta V_{\mathrm{SD}}=2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|} or a \Delta V_{\mathrm{G}}=\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}.

V.4.2 E_{0,S_{0},0}<\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}<E_{0,T^{m}_{1},0}

In this case, we focus on the diamond D and study the current-voltage characteristics for gate voltages $V_{\mathrm{G}}\in[V^{\mathrm{S0}}_{\mathrm{G}},V^{\mathrm{T1}}_{\mathrm{G}}]$ and source-drain bias voltages V_{\mathrm{SD}}\in[0,2\mbox{\displaystyle\left(\tilde{E}^{\alpha}{-1,D^{\sigma}{0},0}\mbox{$\displaystyle\left(V_{\mathrm{G}},0\right)$}-E_{0,S_{0},0}\right)}/\mbox{\displaystyle\left|e\right|}]. At zero e-p coupling, the field-off and field-on conditions result in the same blockade of electron transport but differs in the origin of the blockade. The diamond D is associated with $\displaystyle|\,0,S_{0}\,\rangle$ when optical field is off, but it is related to a partial charged state when the radiation is onFu et al. (2018). The blockade under optical excitation is referred to as anomalous Coulomb blockade in Ref. [​​68].

In contrast with the situation in section V.4.1 where the e-p coupling suppresses the current, in this case, the e-p coupling lifts the anomalous Coulomb blockade, resulting in the photoconductivity within the diamond D. The IV curves in Figure 6b are computed for a representative energy level alignment \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,S_{0},0}+\hbar\omega which is illustrated in Figure 7b. In the bias voltage range [0,2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}], i.e., within the diamond D, Figure 6b shows a current step under optical excitation when $\lambda>0$. The photoinduced current within V_{\mathrm{SD}}\in[0,2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}] does not monotonically change with $\lambda$, but first increases and then decreases. This dependence of photoinduced current on $\lambda$ can be illustrated quantitatively by the analytical solutions.

In order to derive the analytical solutions, we also choose the effective electronic states $\displaystyle|\,-1,D^{\uparrow}\,\rangle$, $\displaystyle|\,-1,D^{\downarrow}\,\rangle$, $\displaystyle|\,0,S_{0}\,\rangle$ and $\displaystyle|\,0,S_{1}\,\rangle$. According to the energy level alignment scheme in Figure 7b, the rate matrix W is constructed as

[TABLE]

The stationary solution then gives the populations as

[TABLE]

and the electric current within the bias voltage range [0,2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}] as

[TABLE]

Similar to the previous section, the zero e-p coupling limit of Eq. (113) reduces to our analytical results in Ref. [​​68], i.e., $I_{\mathrm{S}}=-I_{\mathrm{D}}\rightarrow 0$ and the net charge on the molecule \mbox{\displaystyle\left|Q\right|}=2\mbox{\displaystyle\left|e\right|}P_{-1,D^{\sigma}_{0}}\rightarrow\frac{\mbox{\displaystyle\left|e\right|}k^{\mathrm{field}}}{3k^{\mathrm{field}}+\Gamma} .

Eq. (113) reveals that the photoinduced current can originate from the asymmetry in $\Lambda^{\mathrm{S}}_{0/1}$ and $\Lambda^{\mathrm{D}}_{0/1}$ in the case of E_{0,S_{0},0}<\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}<E_{0,T^{m}_{1},0}. Under the situation specified in Figure 7b, this asymmetry is immediately achieved when a finite $V_{\mathrm{SD}}$ is applied.

It is worthy to point out that the asymmetry in Eq. (113) resembles the Eq. (16) of Ref. [​​62] which attributes the observation of photocurrent to the asymmetric molecule-lead coupling, i.e., $\Gamma_{\mathrm{S}}\neq\Gamma_{\mathrm{D}}$. The Eq. (16) of Ref. [​​62] is derived under an off-resonant tunnelling situation using a HOMO-LUMO model in the single particle picture. The situation of off-resonant tunneling is similar to the energy level alignment considered in this section, especially if we consider that $\Lambda^{\mathrm{\alpha}}_{0}$ and $\Lambda^{\mathrm{\alpha}}_{1}$ are respectively the analogs of the coupling between HOMO and electrode $\alpha$, and the coupling between LUMO and electrode $\alpha$. However, it should be clarified that the asymmetry in Eq. (113) results from the e-p coupling and does not rely on the symmetry of the molecular system, while the photocurrent predicted in Ref. [​​62] requires an intrinsic asymmetry of the molecule.

For the transport characteristics in the case V_{\mathrm{SD}}>2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|} (out of the diamond D), as shown in Figure 6, the IV curves for $\lambda=2.0$ exhibit a transition from the photoconductive behavior to the photoresistive behavior at the fourth current step, in contrast with the IV curves for $\lambda=0.5$ and $1.0$ that exhibit the same photoresistive behavior as in section V.4.1. The current-voltage characteristics for $\lambda=2.0$ thus deserves further discussion. When $\lambda=2.0$, the asymmetry between $\Lambda^{\mathrm{S}}_{0/1}$ and $\Lambda^{\mathrm{D}}_{0/1}$ dominates the transport characteristics within V_{\mathrm{SD}}\in[2,6]\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|} and results in the photoconductive behavior within this bias range. However, as the $V_{\mathrm{SD}}$ increases, the asymmetry gradually disappears due to the activation of more transport channels, which leads to the photoresistive behavior when V_{\mathrm{SD}}>6\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}.

Apart from the representative situation given in Figure 7b, other situations corresponding to \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)})=E_{0,S_{0}}+n\hbar\omega_{\mathrm{vib}} ($n>1$) also deserve further exploration. We highlight these situations because their zero onset voltage of the photocurrent is of experimental interest. The current-voltage characteristics under these situations share a number of key features with the one given in Figure 7b, i.e., (1) e-p coupling induced photocurrent when V_{\mathrm{SD}}\in[0,2n\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}], (2) a photoresistive current-voltage characteristics for $\lambda\leq 1$ and a crossover from photoconductive behavior to photoresistive behavior for $\lambda=2.0$ when V_{\mathrm{SD}}>2n\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}.

In order to gain more insight about the photocurrent under the above situations, we investigate the first current step with the help of Eq. (113). Note that Eq. (113) encodes the dependence of $V_{\mathrm{G}}$ and $V_{\mathrm{SD}}$ in $\Lambda^{\alpha}_{0/1}$ and generally applies to $V_{\mathrm{G}}\in[V^{\mathrm{S0}}_{\mathrm{G}},V^{\mathrm{T1}}_{\mathrm{G}}]$ and V_{\mathrm{SD}}\in[0,2\mathrm{min}(\mbox{\displaystyle\left|V_{\mathrm{G}}-V^{\mathrm{S0}}{\mathrm{G}}\right|},\mbox{\displaystyle\left|V{\mathrm{G}}-V^{\mathrm{T1}}{\mathrm{G}}\right|})]. If we only focus on the gate voltages close to $V^{\mathrm{S0}}_{\mathrm{G}}$ such that (E_{0,S_{1},0}-\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V{\mathrm{G}},V_{\mathrm{SD}}\right)})\gg\hbar\omega_{\mathrm{vib}} and $\Lambda^{\alpha}_{1}\rightarrow 1$, Eq. (113) is simplified to

[TABLE]

Since we are interested in the first current step, i.e., V_{\mathrm{SD}}\in[0,2\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}], Eq. (114) turns to

[TABLE]

where $n$, determined by $V_{\mathrm{G}}$, labels the vibrational excitation associated with $\displaystyle|\,0,S_{0}\,\rangle$ that aligns with \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}).

According to Eq. (115), the photocurrent can be analyzed by the dependence of \mbox{\displaystyle\left|M_{0n}\mbox{$\displaystyle\left(\lambda\right)$}\right|}^{2} on $\lambda$ and $n$, respectively. On one hand, when $n$ is fixed and $\lambda$ is varied, i.e., we are dealing with a specific energy level alignment \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)})=E_{0,S_{0},0}+n\hbar\omega_{\mathrm{vib}}, the photoinduced current approximately maximizes at $\lambda=\sqrt{n}$, which is consistent with the behavior of the first current step in Figure 6b. Figure 8a presents the photocurrent as a function of $\lambda$ for situations of $n=1$ and $n=2$. In Figure 8a, the approximate analytical solution in Eq. (115) agrees well with the numerical calculation and the exact analytical solution in Eq. (113), where a slight deviation is observed under the situation of $n=2$ when $\lambda$ exceeds $2.0$. It should be pointed out that the peak of the $I-\lambda$ curve in Figure 8a slightly deviates from $\sqrt{n}$ due to the monotonous decline of the denominator of Eq. (114) with respect to $\lambda$. On the other hand, when $\lambda$ is fixed and $n$ is varied, the photoinduced current peaks at $n\sim\lambda^{2}$. This explains why the most pronounced conductance lines within the diamond D in Figure 3 shift towards larger gate voltage as the electron phonon coupling increases.

V.4.3 \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,T^{m}_{0},0}

Under this situation, we restrict ourselves to the bias range V_{\mathrm{SD}}\in[0,2\mbox{\displaystyle\left(E_{0,T^{m}{1},0}-E{0,S_{0},0}\right)}/\mbox{\displaystyle\left|e\right|}]. In this bias range, the field-on transport characteristics is activated when $V_{\mathrm{SD}}>0$, while the field-off transport characteristics in this range is completely blocked. The photoinduced current is attributed to the energy level alignment between \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)} and $E_{0,T^{m}_{0},0}$.

The influence of the e-p coupling on the photoinduced transport characteristics is presented in Figure 6c. Among all the explored values of $\lambda$, with the IV curves at $\lambda=0$ as a reference, the e-p coupling suppresses the first three photoinduced currents steps but enhances the subsequent two current steps. These observations are also quantitatively explained by the analytical solution.

Next, we derive the analytical current expression using the effective electronic states $\displaystyle|\,-1,D^{\uparrow}\,\rangle$, $\displaystyle|\,-1,D^{\downarrow}\,\rangle$, $\displaystyle|\,0,S_{0}\,\rangle$, $\displaystyle|\,0,S_{1}\,\rangle$, $\displaystyle|\,0,T^{+1}_{0}\,\rangle$, $\displaystyle|\,0,T^{0}_{0}\,\rangle$ and $\displaystyle|\,0,T^{-1}_{0}\,\rangle$. The vector form of the occupation probabilities is

[TABLE]

According to the energy level alignment illustrated in Figure 7c, the rate matrix W is constructed as

[TABLE]

where we further abbreviate $\Lambda^{\alpha}_{-1,D^{\sigma}_{0}\leftarrow 0,T^{m}_{1}}$ as $\Lambda^{\alpha}_{T}$ for simplicity.

The stationary current corresponding to the situation in Figure 7c is then solved as

[TABLE]

whose zero e-p coupling limit also agrees with the results in Ref. [​​68], i.e.,

[TABLE]

Since the energy level alignment diagram in Figure 7c indicates $\Lambda^{\mathrm{S}}_{T}=\Lambda^{\mathrm{D}}_{T}=0.5\Lambda_{T}$, we can further simplify Eq. (125) into

[TABLE]

with

[TABLE]

and

[TABLE]

where $I_{\mathrm{T}}$ and $I_{\mathrm{as}}$ reveal two origins of the photoinduced current under this situation. First, the component $I_{\mathrm{T}}$ originates from the energy level alignment between \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)} and the triplet states, so it exactly reduces to Eq. (126) as $\lambda\rightarrow 0$ and also contributes to the suppression of the first three photocurrent steps in Figure 6c. Second, the component $I_{\mathrm{as}}$ originates from the asymmetry in $\Lambda^{\mathrm{S/D}}_{0/1}$. $I_{\mathrm{as}}$ resembles Eq. (113) which describes the photocurrent in section V.4.2 and vanishes at $\lambda=0$. It is worth noting that $I_{\mathrm{as}}$ contributes to the enhancement of 4th and 5th photoinduced current steps in Figure 6c.

Figure 8b offers more insights about $I_{\mathrm{T}}$ and $I_{\mathrm{as}}$, where the first photoinduced current step in this case, together with first field-on current step in section V.4.1 is plotted as a function of the e-p coupling $\lambda$. In Figure 8b, the numerical calculations agrees well with the analytical solutions in Eq. (99) and (127). The $I-\lambda$ calculated at $V_{\mathrm{SD}}=0.2$ V and $V_{\mathrm{G}}=V^{\mathrm{T1}}_{\mathrm{G}}$ first displays an exponential decay due to $I_{\mathrm{T}}$ and then exhibits a maximum led by $I_{\mathrm{as}}$ nearby $\lambda=2.0$, while on the contrary the situation of $V_{\mathrm{G}}=V^{\mathrm{S0}}_{\mathrm{G}}$ shows a monotonous decrease in the field-on current with respect to $\lambda$.

Furthermore, the decomposition of the photocurrent into $I_{\mathrm{T}}$ and $I_{\mathrm{as}}$ also explains the role of the e-p coupling on the photoinduced current as shown in Figure 6c. Within V_{\mathrm{SD}}\in[0,6\hbar\omega_{\mathrm{vib}}/\mbox{\displaystyle\left|e\right|}], i.e., during the first three current steps, $I_{\mathrm{T}}$ dominates the photoinduced transport characteristics. In this bias range, $I_{\mathrm{T}}$ is responsible for the e-p coupling induced suppression in the photocurrent, while $I_{\mathrm{as}}$ offers negligible contribution because the asymmetry between $\Lambda^{\mathrm{S}}_{0/1}$ and $\Lambda^{\mathrm{D}}_{0/1}$ is insignificant. The energy level alignment in Figure 7c clearly reveals that the extent of the asymmetry required by $I_{\mathrm{as}}$ grows with $V_{\mathrm{SD}}$. Meanwhile, as $V_{\mathrm{SD}}$ increases, the e-p coupling induced photocurrent suppression due to $I_{\mathrm{T}}$ is also lifted and will be completely eliminated when more transport channels are activated. Henceforth, the gradual elimination of the photocurrent current suppression given by $I_{\mathrm{T}}$, together with the increase in $I_{\mathrm{as}}$, finally results in the e-p coupling enhanced photoinduced current at the 4th and 5th current step in Figure 6c.

V.4.4 \tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}>E_{0,T^{m}_{0},0}

Our discussion for this situation focuses on the bias range V_{\mathrm{SD}}\in[0,2(\tilde{E}^{\alpha}_{-1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}-E_{0,T^{m}_{0},0})/\mbox{\displaystyle\left|e\right|}], beyond which the transport characteristics is governed by Eq. (125). The ranges of $V_{\mathrm{G}}$ and $V_{\mathrm{SD}}$ we specified correspond to the diamond E, where the electron tranport is blocked under both field-off and field-on conditions.

Unlike the photoinduced current described in section V.4.2, the e-p coupling does not lift the blockade under this situation, which has been clearly revealed in the charge stability diagrams. According to the representative energy level alignment diagram in Figure 7d, the transport characteristics in this case is completely blocked by the triplet states. Figure 7d shows that there is no outgoing electronic transitions from triplet states, therefore the triplet states accumulate all the population of the molecular system in the steady state. The analytical solution under this situation is straightforward and obvious, i.e., $P_{0,T^{+1}_{1}}=P_{0,T^{0}_{1}}=P_{0,T^{-1}_{1}}=1/3$. Note that the relevant singlet-triplet transitions is insignificant and negligible because the optical excitations dominate in this situation.

VI Conclusion and Perspective

In conclusion, we have examined the influence of e-p coupling and vibrational relaxation on electron transport through an irradiated molecular junction. Moreover, we have shown that the roles of the triplet states and the energy level alignment between charge states and neutral states are crucial for the transport characteristics. Our main findings are summarized as follows:

(1) In the charge stability diagram, the width of the middle diamond (corresponding to the triplet states) is robust to the e-p coupling. In contrast, the two side diamonds (corresponding to the partial charged states) gradually shrink with the increasing e-p coupling due to the activation of transport channels induced by vibrational transitions. In addition, the effect of the e-p coupling is more pronounced in the unequilibrated phonon regime than in the equilibrated phonon regime because of the vibrational transitions caused by the nonequilibrium vibratinal population.

(2) We have demonstrated a new type of photoconduction mechanism due to the e-p coupling. Our previous study shows that, in the absence of the e-p coupling, the energy level alignment between the charge states and the triplet states uniquely leads to the photocurrent. In this study, in the presence of the e-p coupling, the asymmetry in the rates of charge transfer transitions also results in the photocurrent when the renormalized state energy of a charge state falls between $E_{0,S_{0},0}$ and $E_{0,T^{m}_{1},0}$. This photoconduction mechanism is also significant when the triplet states come into play.

(3) Our study can be used as a guide for the selection of photoconductive molecules. When a high-frequency vibrational mode is coupled to the charge transfer transition, the observation of photocurrent requires $\tilde{E}^{\alpha}_{\pm 1,D^{\sigma}_{0},0}\sim E_{0,T^{m}_{1},0},E_{0,S_{0},n}(n\geq 1)$, i.e., the renormalized state energy of charged states approximately aligns with the triplet states $\displaystyle|\,0,T^{m}_{1}\,\rangle$ or the excited vibrational states of $\displaystyle|\,0,S_{0}\,\rangle$. When a low-frequency vibrational mode dominates, the observation of photocurrent requires $\tilde{E}^{\alpha}_{\pm 1,D^{\sigma}_{0},0}\gtrsim E_{0,S_{0},0}$ or $\tilde{E}^{\alpha}_{\pm 1,D^{\sigma}_{0},0}\lesssim E_{0,T^{m}_{1},0}$, whereas the required proximity of the alignment is determined by the e-p coupling.

(4) Compared with the photoinduced current in the zero e-p coupling limit, the presence of the e-p coupling can either enhance or suppress the photocurrent. In the case of \tilde{E}^{\alpha}_{\pm 1,D^{\sigma}_{0},0}\mbox{\displaystyle\left(V_{\mathrm{G}},0\right)}=E_{0,T^{m}_{1},0}, the influence of the e-p coupling on the photocurrent changes from suppression to enhancement as the source-drain bias increases.

Although we have conducted a comprehensive study on the effect of e-p couplings and vibrational relaxation on an irradiated molecular junction, several issues deserve further discussion. First, the master equation approach cannot include the level broadening due to the molecule-lead coupling and thus fails to describe the tunneling current within the blockade region. Our future work will address this issue through two routes, i.e., the master equation approaches that go beyond 2nd order expansionKoller et al. (2010) and the Hubbard operator NEGF approachGalperin, Nitzan, and Ratner (2008) that intrinsically considers the level broadening at the lowest order expansion. Second, the energy level alignment of molecular systems can be adjusted not only by a gate electrode but also by the functionalization of molecules (e.g., changing the substituent groupLi et al. (2013)), electrochemical gatingHuang et al. (2019) and molecular orbital gatingHines et al. (2010); Choi et al. (2010). Finally, we hope that our work can motivate further studies on irradiated molecular junctions and promote the development of molecular electronics.

Acknowledgements.

This research was supported by Academia Sinica and the Ministry of Science and Technology of Taiwan (MOST 106-2113-M-001-036-MY3). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.

Appendix A The Evaluation of V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b}

For the purpose of evaluating the scattering amplitude V^{*}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b} (V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,NaN-1b}), we recast the molecule-lead coupling Hamiltonian in Eq. (II.1) in terms of single particle basis as,

[TABLE]

where $\hat{d}_{n\sigma}$($\hat{d}_{n\sigma}^{\dagger}$) is the annihilation(creation) operator of an electron with spin $\sigma$ on the single particle level $\displaystyle|\,n\,\rangle$, and t^{*}_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma} (t_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma}) is the scattering amplitude that describes the hopping of an electron with momentum k and spin $\sigma$ from $\displaystyle|\,n\,\rangle$ (electrode $\alpha$) to the electrode $\alpha$ ($\displaystyle|\,n\,\rangle$). In principle, the single particle basis could be any complete orthonormal basis, hereafter we choose the molecular orbitals for convenience.

We next denote a many-body electronic state $\displaystyle|\,N,a\,\rangle$ as $\displaystyle|\,A\,\rangle$ and insert the relationship \hat{1}=\sum_{A}\mbox{\displaystyle|,A,\rangle}\mbox{\displaystyle\langle,A,|} into Eq. (130), which results in

[TABLE]

By comparing Eq. (131) with Eq. (II.1), we obtain the identities

[TABLE]

which relates the V_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,AB}(V^{*}_{\alpha\mbox{\displaystyle\textbf{k}}\sigma,AB}) in the picture of molecular many-electron states with the t_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma}(t^{*}_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma}) in the single particle basis.

Similar to the factorization in Eq. (24), we assume that t_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma}(t^{*}_{n\alpha\mbox{\displaystyle\textbf{k}}\sigma}) could factorize asPeskin (2017)

[TABLE]

where $\zeta_{n\alpha}$($\zeta^{*}_{n\alpha}$) characterizes the overlap between the molecular orbital $\displaystyle|\,n\,\rangle$ and the electrons in the electrode $\alpha$.

Plugging Eq. (134) and (135) into Eq. (132) and (133), we arrive at the definitions of $T_{\alpha,AB}$ and $T^{*}_{\alpha,AB}$ as

[TABLE]

Therefore, it is obvious that $T_{\alpha,AB}$ and $T^{*}_{\alpha,AB}$ survive only when $\displaystyle|\,A\,\rangle$ differs from $\displaystyle|\,B\,\rangle$ by one electron.

In order to explicitly evaluate $T_{\alpha,AB}$ and $T^{*}_{\alpha,AB}$, we approximate the many-electron states considered in the main text in terms of the Kohn-Sham orbitals,

[TABLE]

where $\hat{c}^{\dagger}_{n\sigma}$ and $\hat{c}_{n\sigma}$ are respectively the creation operator and the annihilation operator of an electron on Kohn-Sham orbital $\displaystyle|\,n\,\rangle$ with spin $\sigma$ , H refers to HOMO, L refers to LUMO, and $\displaystyle|\,\chi\,\rangle$ is a reference state with empty HOMO and LUMO. Here we assume that all the many-electron states could be constructed from the same set of Kohn-Sham orbitals.

Note that, in $T_{\alpha,N\pm 1b,Na}$, only up to one term survives in the summation over orbital index $n$, i.e. either $\zeta_{1\alpha}$ or $\zeta_{2\alpha}$. We suppose that $\zeta_{1\alpha}=\zeta_{2\alpha}=\zeta_{\alpha}$, and then reorganize \frac{2\pi}{\hbar^{2}}\mbox{\displaystyle\left|T_{\alpha,N\pm 1b,Na}\right|}^{2}\bar{J}_{\alpha} into $\Gamma_{\alpha}\nu_{N\pm 1b,Na}$, where the characteristic rate of charge transfer transition $\Gamma_{\alpha}$ is defined as

[TABLE]

and the dimensionless coupling coefficient $\nu_{N\pm 1b,Na}$ of the transition \mbox{\displaystyle|,N\pm 1,b,\rangle}\leftrightarrow\mbox{\displaystyle|,N,a,\rangle} has been given in Table S1 of Ref. [​​68].

Appendix B Derivation of $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ in Eq. (III.4)

According to Eq. (III) and (III.4), $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ is defined as

[TABLE]

where \rho_{0S_{0}\nu_{0},0S_{1}\nu_{1}}=\mbox{\displaystyle\langle,0,S_{0},\nu_{0},|}\hat{\rho}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,0,S_{1},\nu_{1},\rangle}.

In order to derive $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$, we follow the same procedure as we used to derive the rate equation. Since the SOC-induced transitions are significantly slower than the other processes, we restrict ourselves to contributions from $\hat{H}_{\mathrm{m-l}}$ and $\hat{H}_{\mathrm{m-th}}$, which are hereafter referred to as $\kappa^{\mathrm{m-l}}_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ and $\kappa^{\mathrm{m-th}}_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ respectively, i.e.,

[TABLE]

The derivation of $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ starts from taking the off-diagonal matrix element of Eq. (III.2). Unlike the derivation of charge transfer rate equation, not all the terms in Eq. (III.2) contribute to $\kappa^{\mathrm{m-th}}_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$. In the following, we take the derivation for the 1st and 4th terms in Eq. (III.2) as examples, one of which contributes to $\kappa_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$ and the other does not.

The matrix element of the $1$st term in Eq. (III.2) with respect to $\displaystyle|\,N,a,\nu_{a}\,\rangle$ and $\displaystyle|\,M,b,\nu_{b}\,\rangle$ is derived as

[TABLE]

Similar to the Section III.2, the secular approximation leads to

[TABLE]

which means that Eq. (B) does not contribute to $\kappa^{\mathrm{m-l}}_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$. Note that we here neglect the possibility of $E_{N,a,\nu_{a}}-E_{M,b,\nu_{b}}=E_{M-1,b_{2},\nu_{2}}-E_{N-1,b_{1},\nu_{1}}$ when $E_{N,a,\nu_{a}}\neq E_{M,b,\nu_{b}}$, which is a reasonable assumption when considering the electronic structure of a molecular system.

We next derive the matrix element of the $4$th term in Eq. (III.2), which gives

[TABLE]

With the help of the secular approximation, we have \mbox{\displaystyle\mathbf{e}^{\frac{i}{\hbar}\mbox{$\displaystyle\left(E_{N,a,\nu_{a}}-E_{N,b_{1},\nu_{1}}\right)$}t}}\rightarrow\delta_{a,b_{1}}\delta_{\nu_{a},\nu_{1}} , which leads to

[TABLE]

where $\kappa^{\mathrm{m-l}(4)}_{Na\nu_{a},Mb\nu_{b}}$ refers to the contribution from the 4th term in Eq. (III.2). In order to find out the physical meaning of $\kappa^{\mathrm{m-l}(4)}_{Na\nu_{a},Mb\nu_{b}}$, we apply the wide band limit and the formula

[TABLE]

where $\cal P$ refers to the Cauchy principal value. As a result, we obtain

[TABLE]

where

[TABLE]

is an energy shift caused by the charge transfer transition \mbox{\displaystyle|,N,a,\nu_{a},\rangle}\rightarrow\mbox{\displaystyle|,N-1,b_{2},\nu_{2},\rangle}.

We can similarly evaluate the other terms in Eq. (III.2) and then obtain

[TABLE]

Next, we move forward onto the derivation of $\kappa^{\mathrm{m-th}}_{0S_{0}\nu_{0},0S_{1}\nu_{1}}$, which starts from taking the off-diagonal matrix element of Eq. (III.3). We then obtain

[TABLE]

where the terms \mbox{\displaystyle\langle,N,a,\nu_{a}-1,|}\hat{\rho}^{I}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,M,b,\nu_{b}-1,\rangle} and \mbox{\displaystyle\langle,N,a,\nu_{a}+1,|}\hat{\rho}^{I}\mbox{\displaystyle\left(t\right)}\mbox{\displaystyle|,M,b,\nu_{b}+1,\rangle} are dropped because they introduce higher order terms of the coupling Hamiltonian $\hat{H}_{\mathrm{m-th}}$, which is reasonable in the limit of weak system-bath coupling. According to Eq. (III.3), we can express $\kappa^{\mathrm{m-th}}_{0S_{0}\nu_{a},0S_{1}\nu_{b}}$ as

[TABLE]