Characterization and control of matter by optical means is at the forefront of research both due to fundamental insights and technological promise. Theoretical modeling of periodically driven systems is a prerequisite to understanding and engineering nanoscale quantum devices for quantum technologies. Here, we develop a theory for transport and optical response of molecular junctions, open nonequilibrium quantum systems, under external periodic driving. Periodic driving is described using the Floquet theory combined with nonequilibrium Green’s function description of the system. Light–matter interaction is modeled by employing the self-consistent Born approximation. A generic three-level model is utilized to illustrate the effect of the driving on optical and transport properties of junctions.

Recent progress in laser techniques combined with technological advances at the nanoscale resulted in the possibility of optical measurements in current carrying single-molecule junctions. In particular, bias induced fluorescence1–3 and phosphorescence4 were reported in the literature. Electroluminescence allowed to visualize intra-molecular interactions,5,6 explore real space energy transfer,7 and probe charge fluctuations in biased molecular junctions.8 Raman spectroscopy yields information on the vibrational structure and electron flux induced heating of vibrational and electronic degrees of freedom in single-molecule junctions.9–12 Recently, experimental measurements of the ultra-strong light–matter interaction in single molecule cavities (so far without current) were reported in the literature.13,14 The combination of molecular electronics and optical spectroscopy evolved into a new direction of research coined molecular optoelectronics.15,16

Control of response in nanoscale systems by external driving is an active area of research both experimentally and theoretically due to technological promise of engineering and control in quantum devices. For example, periodic radio frequency potential modulations in a suspended carbon nanotube junction facilitated experimental studies of strong coupling between tunneling electron and nanomechanical motion.17 AC-driven charging and discharging of quantum dots were employed in experimental verification of the existence of quantum stochastic resonance, where intrinsic fluctuations lead to amplification and optimization of a weak signal.18 In single molecule junctions, periodic driving by a set of laser pulse pairs was suggested as a tool to study sub-picosecond intra-molecular dynamics.19,20

Here, we consider a molecular junction driven by a time-periodic external field. The effective electronic properties of this driven system are characterized by developing a theory of its optical and transport response in this highly nonequilibrium setting. Theoretical studies of periodic driving in junctions usually employ the quantum master equation21–25 and nonequilibrium Green’s function26–28 approaches. Periodicity in external driving allows mapping of the original time-dependent problem into effective time-independent formulation in an extended Floquet space.29 The Floquet theory was used in numerous theoretical studies of external driving on transport in junctions. Studies employing a combination of the theory with the Schrödinger equation,30 scattering matrix approach,31 quantum master equation,32–35 and Green’s functions36–40 are available in the literature. For strongly correlated systems, combinations of the Floquet with the slave boson,41 dynamical mean field theory,42 and functional renormalization group43 were also formulated. Reference 44 considered the periodically driven strongly interacting Fermi–Hubbard model employing the generalization of the Schrieffer–Wolf transformation to periodically driven systems using Floquet theory.

Note that for molecular junctions, the characteristic strength of molecular coupling to contacts Γ ∼ 0.01–0.1 eV45 is of the same order of magnitude or stronger than thermal energy kBT ∼ 0.01 eV. Thus, the Redfield/Lindblad quantum master equation, which relies on assumption Γ ≪ kBT, is not suitable and the utilization of the nonequilibrium Green’s function (NEGF) technique is preferable. We consider a case where the light–matter interaction is smaller than the system-bath coupling so that diagrammatic perturbation theory can be employed and treat the interaction within the self-consistent Born approximation (SCBA). Periodicity of driving makes the Floquet theory useful.

Here, we formulate a Floquet–NEGF–SCBA theory of transport and optical response of single-molecule junctions under external periodic driving. The formulation complements the study of optical properties of current carrying junctions46 by accounting for periodic external driving and extends the recent study on optical absorption properties of driven matter47 to the realm of open nonequilibrium molecular systems. We note in passing that the second order of the diagrammatic expansion employed in the SCBA makes our consideration somewhat similar to that of Ref. 44 because results of the Schrieffer–Wolf transformation can be effectively achieved within second order perturbation theory.48 

The structure of the paper is as follows: In Sec. II, we introduce the junction model and formulate the Floquet version of the NEGF–SCBA treatment of electron and photon fluxes in the junction. We use the Floquet–NEGF–SCBA to illustrate effects of quantum coherence on junction responses in Sec. III. Section IV concludes and outlines directions for future research.

Here, we present a model of a junction and briefly describe the nonequilibrium Green’s function (NEGF) approach for the simulation of photon and electron fluxes, where we treat the light–matter interaction at the level of the self-consistent Born approximation (SCBA). After this, we present the Floquet formulation of the NEGF–SCBA problem.

We consider a junction consisting of a molecule M coupled to two contacts L and R and to quantum modes {α} of the radiation field pt. For simplicity, the molecule is modeled as a non-interacting (mean-field or DFT) system. Note that this is the usual49–51 although not always a reliable52 level of modeling in first principles simulations. In addition, the molecule is subjected to external driving by monochromatic classical field E(t),

E(t)=E0cos(ω0t+ϕ0).
(1)

Contacts are modeled as reservoirs of free electrons each at its own equilibrium (see Fig. 1). The external driving distorts the electronic structure of the junction leading to effective properties that can, in principle, be very different from those observed at equilibrium.

FIG. 1.

Sketch of a molecular junction subjected to external driving.

FIG. 1.

Sketch of a molecular junction subjected to external driving.

Close modal

The Hamiltonian of the model is (here and below = e = kB = 1)

Ĥ(t)=ĤM(t)+B=L,R,ptĤB+V^MB,
(2)

where ĤM and ĤB (B = L, R, pt) represent the decoupled molecule and baths (contacts and radiation field), while V^MB are the couplings. Explicit expressions are

ĤM(t)=m1,m2MHm1m2Mμm1m2E(t)d^m1d^m2,ĤL(R)=kL(R)εkĉkĉk;  Ĥpt=αωαâαâα,V^ML(R)=mMkL(R)Vkmĉkd^m+H.c..V^M,pt=m1,m2MαVα,m1m2ptâαd^m1d^m2+H.c..
(3)

Here, d^m (d^m) and ĉk (ĉk) create (annihilate) the electron in orbital m of the molecule and state k of the contacts, respectively. âα (âα) creates (destroys) excitation quanta in mode α of the radiation field. HM is the molecular (or Kohn–Sham) orbitals basis representation of the effective single particle Hamiltonian of the molecule. εk and ωα are the energy of the free electron in state k in contacts and the frequency of the radiation field mode α, respectively. Vkm is the matrix element for electron transfer from molecular orbital m to state k in contacts. Vα,m1m2pt is the matrix element for intra-molecular (optical) electron transfer from orbital m2 to orbital m1 while creating excitation quanta (photon) in mode α. We stress that coupling to the radiation field in (3) is completely general, i.e., it is not restricted to the rotating wave approximation.

Below, for simplicity, we assume

Vm1m2,αpt=Um1m2[Vαpt]*,
(4)

where U is the matrix in molecular orbital space with 1 for allowed optical transitions and 0 otherwise.

The nonequilibrium Green’s function (NEGF) in the subspace of molecular orbitals is defined on the Keldysh contour as

Gm1m2(τ1,τ2)=iTcd^m1(τ1)d^m2(τ2),
(5)

where Tc is the contour ordering operator and τ1,2 are the contour variables. Retarded (r) and lesser (<)/greater (>) projections of the Greens’ function (5) are usually obtained by solving differential (or Kadanoff–Baym) and integral (or Keldysh) forms of the Dyson equation, respectively,53,54

iIt1HMμE(t1)Gr(t1,t2)  dtΣr(t1,t)Gr(t,t2)=Iδ(t1t2),
(6)
G(t1,t2)=dt3dt4Gr(t1,t3)Σ(t3,t4)Ga(t4,t2).

Here, Hamiltonian HM, molecular dipole moment μ, Green’s function G, and self-energy Σ projections are written as matrices in the subspace of molecular orbitals. I is the unity matrix in the subspace and Ga(t4,t2)[Gr(t2,t4)] is the advanced projection of the Green’s function (5). The self-energy accounts for contributions of the baths (contacts and radiation field) into molecular dynamics,

Σm1m2(τ1,τ2)=B=L,R,radΣm1m2B(τ1,τ2),Σm1m2L(R)(τ1,τ2)=kL(R)Vm1kgk(τ1,τ2)Vkm2,Σm1m2pt(τ1,τ2)=m3,m4MGm3m4(τ1,τ2)Um3m1F(τ1,τ2)Um4m2+Um4m2F(τ2,τ1)Um1m3,
(7)

where

F(τ1,τ2)α|Vαpt|2fα(τ1,τ2).
(8)

gk(τ1,τ2)iTcĉk(τ1)ĉk(τ2)0 is the Green’s function of the free electron in state k, and fα(τ1,τ2)iTcâα(τ1)âα(τ2)0 is the Green’s function of the photon in mode α of the radiation field.

Note that while expressions for self-energies due to contacts L and R are exact, because of non-quadratic molecule–radiation field coupling in (3), the self-energy due to the radiation field can be treated only approximately. We utilize second order diagrammatic expansion (the self-consistent Born approximation) for Σpt (for more details, see Ref. 55). Also note that for simplicity, we disregard the effect of the molecule on radiation field modes.

Once Green’s function is known, it can be used for the evaluation of molecular responses to external perturbations. Below, we are interested in simulating photon flux Ipt(t)55 and electron current IL(R)(t) through interface L (R),56 

Ipt(t)=2RetdtF<(t,t)Π>(t,t)F>(t,t)Π<(t,t),IL(R)(t)=2Retdt TrΣL(R)<(t,t)G>(t,t)ΣL(R)>(t,t)G<(t,t).
(9)

Here, the trace is over molecular subspace and Π is the lesser/greater projection of the photon self-energy due to coupling to electrons evaluated within the self-consistent Born approximation,

Π(τ1,τ2)=i TrUG(τ1,τ2)UG(τ2,τ1).
(10)

Note that in (9), we follow the convention where the positive photon flux goes from the molecule to the radiation field, while the positive electron flux goes from contacts into the molecule.

While for general time-dependent drivings, one has to solve the Dyson equation (6) numerically on a two-dimensional time grid,57 periodic driving allows us to formulate an effective time-independent problem in an extended Floquet space.29 Because of time periodicity of the Hamiltonian, one-sided Fourier transform of the Green’s function retarded projection can be expanded in discrete series of Floquet modes

Gr(t1,t2)=f=+dE2πGr(f;E)eiE(t1t2)+ifω0t1.
(11)

Substituting this expansion into the Dyson equations (6) leads to

I(Efω0)HMGr(f;E)+12μE0eiϕ0Gr(f1;E)  +eiϕ0Gr(f+1;E)f1,f2=+Σr(f1,f2;E[f+f1]ω0)  ×Gr(f+f1+f2;E)=Iδf,0,
(12)
G(f1,f2;E)=f3,f4=+Gr(f1f3;Ef3ω0),×Σ(f3,f4;E)Ga(f2f4;Ef4ω0).

Here, Ga(f;E)[Gr(f;E)],

G(t1,t2)=f1,f2=+dE2πG(f1,f2;E)eiE(t1t2)+if1ω0t1if2ω0t2,
(13)

and explicit expressions for the self-energy projections in the Floquet space are given in  Appendix A. Equation (12) presents the effective time-independent formulation of the Floquet-NEGF-SCBA.

Once Floquet space projections of the Green’s function are known, we use them to evaluate dc (period averaged) photon and electron fluxes. Substituting (11) and (13) to (9) yields

Iptdc=f=+dω2πF<(ω)Π>(f,f;ω+fω0)F>(ω)Π<(f,f;ω+fω0),IL(R)dc=f=+dE2π TrΣL(R)<(Efω0)G>(f,f;E)ΣL(R)>(Efω0)G<(f,f;E).
(14)

Here, G<(f, f; E) and G>(f, f; E) are calculated from (12), expressions for Π<(f, f; ω) and Π>(f, f; ω) are given in (A7),

F<(ω)=iγ(ω)Npt(ω),F>(ω)=iγ(ω)Npt(ω)+1,ΣL(R)<(E)=iΓL(R)fL(R)(E),ΣL(R)>(E)=iΓL(R)[1fL(R)(E)],
(15)

where

γ(ω)2πα|Vαpt|2δ(ωωα),Γm1m2L(R)(E)2πkL(R)Vm1kVkm2δ(Eεk)
(16)

are the radiative energy dissipation rate and electronic decay rate, respectively, fL(R)(E) is the Fermi–Dirac distribution, and Npt(ω) is the population of radiation field modes at frequency ω.

We note that optical absorption theory presented in Ref. 47 is the zero bias quasiparticle limit of the present formulation with focus on incoming flux and neglecting self-consistency in the treatment of light–matter interaction (see  Appendix B for details). In the case of quantum treatment of radiation field, the latter assumption is known to violate conservation laws58,59 and may lead to qualitative failures.60 Thus, it is not applicable in optical studies of open quantum systems.61 Below, we present numerical illustrations of the developed Floquet–NEGF–SCBA formulation.

Simulations start from the evaluation of Green’s functions employing (12) in the absence of radiation field (i.e., taking self-energy Σpt to be zero). Resulting expressions are utilized to evaluate self-energy due to coupling to the radiation field Σpt—last two lines in Eq. (7). Thus, the obtained self-energy is used in the simulation of the updated Green’s functions employing (12). This self-consistent procedure continues until convergence. The solution was assumed to reach convergence when the difference of level populations for each Floquet mode, nm(f)idEGmm<(f,f;E)/2π, at subsequent steps, s and s + 1, of the iterative procedure is less than 10−6: |nm(s+1)(f)nm(s)(f)|<106 for each f.

In our simulations, we use a generic three level model with each level independently coupled to the contacts,

Hm1m2M=δm1,m2εm1,Γm1m2L(R)=δm1,m2Γm1L(R).
(17)

For simplicity, we also assume the wide band approximation for contacts and radiation bath, i.e., ΓL(R) and γ are energy independent and corresponding Lamb shifts, ΛL(R) and λ, are zero. Parameters and results of the simulations are presented in terms of the arbitrary unit of energy U0 and the corresponding unit of flux I0 = U0/. Unless stated otherwise, parameters are the following: kBT = 0.03, electron escape rates are ΓmL=ΓmR=0.01, and the energy dissipation rate is γ = 10−4. Population of the modes corresponding to the pumping frequency ωp is N(ωp) = 1, and other modes are assumed to be unpopulated. The phase of the driving field ϕ0 = 0. Fermi energy is taken as origin, EF = 0, and bias Vsd is applied symmetrically μL = EF + |e|Vsd/2 and μR = EF − |e|Vsd/2. The Floquet space in the simulations was restricted to 21 modes, i.e., f ∈{−10, … , 0, … , 10}. Simulations are performed on the energy grid spanning the range from −2 to 2 with step 10−3.

We start from the model considered in Ref. 47. Here, ε0 = −0.5, ε1 = 0.5, and ε2 = 1.1. The radiation field is assumed to be coupled to optical transition between orbitals 0 and 1, while external driving couples orbitals 1 and 2 and is set to be in resonance with the transition, ω0 = 0.6 [see Fig. 2(a)]. As discussed in Ref. 47, at zero applied bias voltage (Vsd = 0), the single peak in absorption spectrum splits into two for stronger driving [compare solid and dashed lines in Fig. 2(b)]. The splitting starts at μE0m=12Γm (here ΓmΓmL+ΓmR) and increases with strength of the driving [see Fig. 2(c)].

FIG. 2.

Three-level junction responses to external drivings. (a) The sketch of the junction, (b) photon flux vs pumping frequency ωp for μE0=0 (blue solid line) and 0.1 (blue dashed line) at Vsd = 0, (c) map of the photon flux vs ωp and μE0 at Vsd = 0, (d) photon flux for μE0=0.1 at Vsd = 0 (blue dashed line) and 1 (blue dotted line), (e) electron density of states for μE0=0 (blue solid line) and 0.2 (red dashed line), and (f) electron flux vs bias Vsd for μE0=0 (blue solid line) and 0.2 (red dashed line). Red circles in (b) and (d) show the results of the Floquet–NEGF–SCBA simulations. Other results are obtained by using the rotational frame of the field (see the text for parameters and further details).

FIG. 2.

Three-level junction responses to external drivings. (a) The sketch of the junction, (b) photon flux vs pumping frequency ωp for μE0=0 (blue solid line) and 0.1 (blue dashed line) at Vsd = 0, (c) map of the photon flux vs ωp and μE0 at Vsd = 0, (d) photon flux for μE0=0.1 at Vsd = 0 (blue dashed line) and 1 (blue dotted line), (e) electron density of states for μE0=0 (blue solid line) and 0.2 (red dashed line), and (f) electron flux vs bias Vsd for μE0=0 (blue solid line) and 0.2 (red dashed line). Red circles in (b) and (d) show the results of the Floquet–NEGF–SCBA simulations. Other results are obtained by using the rotational frame of the field (see the text for parameters and further details).

Close modal

It is easy to understand the effect by using the rotating frame of the driving field. Because driving is at resonance, ω0 = ε2ε1, utilization of the rotating wave approximation (RWA) allows us to map the original time-dependent problem Ĥ(t), Eq. (2), onto an effective Hamiltonian H¯^(t),

H¯^(t)=iddteŜ(t)eŜ(t)+eŜ(t)Ĥ(t)eŜ(t),Ŝ(t)=iω02tn^0+n^1n^2,
(18)

where n^md^md^m. The transformation leads to

H¯^(t)=m=02ε¯md^md^m12μE0d^1d^2eiϕ0+d^2d^1eiϕ0+kεkĉkĉk+αωαâαâα+m=02kV¯km(t)ĉkd^m+H.c.+α[Vαpt]*âαd0d^1+H.c..
(19)

Here, ε¯0,1=ε0,1+ω0/2, ε¯2=ε2ω0/2, V¯k0(1)(t)=Vk0(1)eiω0t/2, and V¯k2(t)=Vk2eiω0t/2. Note that in (19), the external driving enters as a time-independent electron hopping between renormalized degenerate levels 1 and 2, ε¯1=ε¯2. Moreover, because the original levels are well separated, |εm1εm2|Γm1,Γm2, bath induced coherences among them can be safely disregarded. In this regime, the transition to the rotating frame of the field allows mapping onto an effective time-independent problem. Indeed, for diagonal electron dissipation matrix Γ, time factors in V¯km(t) will result in a shift of chemical potentials compensating the renormalization of the orbital levels, i.e., μL(R) is shifted up by ω0/2 for levels 0 and 1 and down by ω0/2 for level 2. This time-independent problem can be treated within the standard NEGF–SCBA.62 

In the rotating frame of the field, the splitting, μ0E0, is a manifestation of absorption into two eigenlevels resulting from the diagonalization of the 1–2 block of the molecular Hamiltonian. The two peaks in the spectrum become distinguishable when the splitting is bigger than the broadening of the eigenlevels. Note that for the parameters of the simulation, the RWA is quite a reasonable approximation. Effects beyond the RWA, such as difference in absorption peak heights [dips in Fig. 2(b)], are small—compare the Floquet–NEGF–SCBA treatment of the original problem (circles) with NEGF–SCBA results for effective time independent model (dashed line). At finite bias, when the depopulation of level 0 and the population of level 1 by electron transfer between the molecule and contacts become possible, absorption is substituted by emission, dips become less pronounced, and depths differ due to a different population of the eigenlevels at Vsd = 1 [compare dashed and dotted lines in Fig. 2(d)].

Similarly, splitting can be observed in transport characteristics of the junction. Figure 2(e) shows splitting in the density of states A (conductance) induced by μE0=0.2 external driving (compare solid and dashed lines), while Fig. 2(f) demonstrates consequences of the splitting on current–voltage characteristics of the junction. The difference is due to tunneling via driven levels of the model. The calculations are performed in the rotating frame of the field.

When bath-induced coherences are non-negligible, consideration becomes more involved. In this case, one cannot formulate an effective time-independent problem, and the Floquet space consideration becomes important. For example, such a situation happens when the radiation field is coupled to two optical transitions with external driving in between levels in the transitions. We consider a three-level model with ε0 = −0.1, ε1 = 0.1, and ε2 = 1. The radiation field is assumed to be coupled to optical transitions between orbitals 0 and 2 and between 1 and 2, while external driving couples orbitals 0 and 1 [see Fig. 3(a)]. As previously, the driving frequency is taken at resonance with its optical transition, ω0 = ε1ε0 = 0.2. Contrary to the model in Fig. 2(a), where the single peak splits into two on increasing the driving strength, absorption in Fig. 3(b) shows transition from one peak to split spectrum (dashed to solid to dotted lines in the panel b).

FIG. 3.

Photon flux vs pumping frequency ωp in a three-level junction model at Vsd = 0. (a) The sketch of the junction, (b) photon flux for several driving strengths μE0 at ω0 = 0.2, and (c) photon flux for several driving frequencies ω0 at μE0=0.2. Simulations are performed with the Floquet–NEGF–SCBA approach (see the text for parameters and further details).

FIG. 3.

Photon flux vs pumping frequency ωp in a three-level junction model at Vsd = 0. (a) The sketch of the junction, (b) photon flux for several driving strengths μE0 at ω0 = 0.2, and (c) photon flux for several driving frequencies ω0 at μE0=0.2. Simulations are performed with the Floquet–NEGF–SCBA approach (see the text for parameters and further details).

Close modal

Because the radiation field couples orbitals 0 and 1, formulation of the effective time-independent problem becomes impossible. Nevertheless, transition to the rotating frame similar to the one above yields a qualitative understanding of the observed effects of quantum coherence on the molecular absorption spectrum. One can show that the absorption spectrum of the model in Fig. 3(a) having a single peak (dip in Fig. 3) at ε2ε0 in the absence of driving when the driving is present can have four peaks at frequencies (see  Appendix C for details)

ω=ε2ε0+ε1±ω02±12(ε1ε0ω0)2+μE02.
(20)

In particular, for the parameters of calculation, one expects three peaks for μE0=ω0 and four peaks otherwise [two additional peaks for each μE0 are outside of the ωp range in Fig. 3(b)]. Heights of the peaks and their broadenings can only be obtained from the simulations. Another experimentally easily accessible parameter is the driving frequency. Here, qualitatively similar behavior is also observed [see Fig. 3(c)].

Finally, we demonstrate the coherent control of absorption and transport in a junction beyond RWA considerations. The molecule is presented by its HOMO ε0 = −0.5, LUMO ε1 = 0.4, and LUMO+1 ε2 = 0.6 levels. External driving couples LUMO and LUMO+1 and is chosen to be in resonance with the optical transition, ω0 = 0.2. Its coupling strength is taken as μE0=0.2. The radiation field couples the ground state (HOMO) with the two excited states (LUMO and LUMO+1). The energy dissipation rate is γ = 0.01. Coherence between radiation and driving fields results in control of the junction responses. In particular, changing the driving field phase ϕ0 can affect optical absorption, which, in turn, affects the current through the junction [compare solid and dashed lines in Fig. 4(b) for absorption and in Fig. 4(c) for current]. Note that the negative sign of photon flux indicates absorption, while the negative sign of electron flux indicates electron transfer from the molecule to the contact R. Also note that possibility of coherent control of optical properties was first discussed in Ref. 63.

FIG. 4.

Coherent control in a three-level junction model at Vsd = 0.5. (a) The sketch of the junction and (b) photon Ipt and (c) electron IR fluxes [Eq. (14)] for ϕ0 = 0 (solid line, triangles, blue) and ϕ0 = π/2 (dashed line, circles, red). Simulations are performed with the Floquet–NEGF–SCBA approach (see the text for parameters and further details).

FIG. 4.

Coherent control in a three-level junction model at Vsd = 0.5. (a) The sketch of the junction and (b) photon Ipt and (c) electron IR fluxes [Eq. (14)] for ϕ0 = 0 (solid line, triangles, blue) and ϕ0 = π/2 (dashed line, circles, red). Simulations are performed with the Floquet–NEGF–SCBA approach (see the text for parameters and further details).

Close modal

Characterization and control of matter by optical means is at the forefront of research both due to fundamental insights and technological promise. For example, in molecular junctions, periodic driving of the system by laser pulse pairs was proposed as a tool to study intra-molecular sub-picosecond dynamics in biased junctions. Theoretical modeling of periodically driven systems is a prerequisite to understanding and engineering nanoscale quantum devices for quantum technologies.

We consider the open nonequilibrium optoelectronic system under external periodic driving. The system was treated within the nonequilibrium Green’s function (NEGF) approach, and the light–matter interaction was assumed to be small enough to allow treatment at the second order of the diagrammatic expansion—the self-consistent Born approximation (SCBA). Periodic driving is accounted for within the Floquet theory, which maps the original time-dependent Hamiltonian onto the effective time-independent problem in the extended Floquet space.

Combined Floquet–NEGF–SCBA expressions were formulated, and the approach was illustrated with generic three-level junction model simulations, demonstrating coherent control of optical response in such devices. The data that support the findings of this study are available from the corresponding author upon reasonable request. The formulation complemented the study of optical properties of current carrying junctions proposed in Ref. 46 by accounting for periodic external driving and extended the recent study on optical absorption properties of matter put forward in Ref. 47 to the realm of open nonequilibrium molecular systems. Incorporation of strong intra-system interactions and application of the methodology in ab initio simulations are goals for future research.

This work was supported by the National Science Foundation under Grant Nos. CHE-1565939 (M.G.) and CHE-1553939 (I.F.). G.C. acknowledges the UCSD Physical Sciences Undergraduate Summer Research Award.

Here, we give expressions for the self-energies utilized in (12). The expressions are obtained by substituting (11) and (13) into corresponding projections of (7). The retarded projection is

Σr(f1,f2;E)=δf1,0δf2,0ΣL(R)r(E)+Σptr(f1,f2;E),
(A1)
ΣL(R)r(E)=ΛL(R)(E)i2ΓL(R)(E),
(A2)
Σptr(f1,f2;E)=idω2πUG>(f1,f2;Eω)Fr(ω)+δf2,0Gr(f1;Eω)F<(ω)U+UG<(f1,f2;E+ω)Fa(ω)+δf2,0Gr(f1,E+ω)F<(ω)U,
(A3)

where F<(ω) is defined in (15),

Fr(ω)=λ(ω)i2γ(ω),Fa(ω)=Fr(ω)*,
(A4)

γ(ω) and ΓL(R)(E) are defined in (16), and λ(ω) and ΛL(R)(E) are the real parts of the corresponding retarded projections associated with γ(ω) and ΓL(R)(E) via the Kramers–Kronig relations.

Lesser and greater projections are

Σ(f1,f2;E)=δf1,0δf2,0ΣL(R)(E)+Σpt(f1,f2;E),
(A5)

where ΣL(R)(E) are defined in (15) and

Σpt(f1,f2;E)=idω2πUG(f1,f2;Eω)F(ω)U+UG(f1,f2;E+ω)F(ω)U.
(A6)

Finally, substituting (13) into lesser and greater projections of (10) leads to

Π(f1,f2;ω)=if3,f4f5,f6=+δf1,f3f5δf2,f4f6dE2π Tr×UGf3,f4;E+ω2UGf6,f5;Eω2.
(A7)

Here, we compare our consideration to the theory of optical absorption presented in Ref. 47.

We start from the expression for the absorption rate given in Eq. (14) of the latter reference. In terms of our definitions, the expression is

I(ωp)=limt|Vppt|2221tt0t0tdt1t0tdt2Reeiωp(t1t2)+eiωp(t1+t2)μ^(t2)μ^(t1),
(B1)

where

μ^=m1,m2MUm1m2d^m1d^m2,
(B2)

Vppt and Um1m2 are defined in Eq. (4) and time-dependence in the correlation function of (B1) is with respect to the Hamiltonian excluding coupling to the radiation field (interaction picture). Note that the latter results in a non-conserving expression and is applicable only in cases when coupling to the radiation field is small (relative to other energy scales of the problem) so that mistake in the treatment is numerically negligible. Note also that a factor of 2 in the denominator in (B1) comes from assumed cosine time dependence of the classical consideration of the radiation field in Ref. 47, which is effectively included in Vppt in our quantum consideration. Thus, it will be dropped below. Similarly, we drop eiωp(t1+t2) because, as was discussed in Ref. 47, the term does not contribute to the absorption spectrum. Finally, we take = 1 in our consideration.

Implementing the adjustments and using the fact that the integrand is a symmetric function of t1 and t2, Eq. (B1) can be rewritten as

I(ωp)=limt|Vppt|21tt0t0tdt1t0t1dt22Reeiωp(t1t2)μ^(t2)μ^(t1).
(B3)

Below, we show that this expressions yields time-averaging of in-scattering (absorption) photon flux, first term in Eq. (9), i.e., we show that

I(ωp)=limt1tt0t0tdt1Iptabs(t1),
(B4)

where

Iptabs(t1)=2Ret0t1dt2Fp<(t1,t2)Π>(t2,t1).
(B5)

Note that in our consideration, t0 was taken at − and the convention of negative sign for the in-coming photon flux was followed. Also note that for the dc part of the flux (the only part considered here and in Ref. 47), averaging (B4) is redundant.

Taking into account that F< in (9) is defined as the sum over modes of the radiation field [see Eq. (8)] to get (B5), one has to consider a single mode of the radiation field with frequency ωp. From definition (15) and assuming singly populated mode Npt(ωp) = 1, Green’s function Fp< for the mode can be written in time domain as

Fp<(t1,t2)=i|Vppt|2eiωp(t1t2).
(B6)

Thus, comparing (B3) to (B5) and (B6), one concludes that

Π>(t2,t1)=iμ^(t2)μ^(t1),
(B7)

The right-hand side of this expression is the definition of the two-particle Green’s function. When coupling to the radiation field is neglected (as is done in Ref. 47) and taking into account that the rest of the Hamiltonian (3) does not contain many-body interactions, one can employ Wick’s theorem expressing two-particle Green’s function in terms of the product of single-particle GFs. In particular, employing Wick’s theorem in the right hand side of (B7) results in an expression similar to greater projection of (10). The crucial difference is that the resulting expression contains zero-order single-particle GFs, while Eq. (10) employs full (dressed in coupling to radiation field) GFs. This completes the derivation of theory presented in Ref. 47 from the Floquet–NEGF–SCBA consideration.

Here, we present the qualitative analysis leading to Eq. (20).

The model shown in Fig. 3(a) has Hamiltonians (2) and (3) with Hm1m2M=δm1,m2εm1 (m ∈{0, 1, 2}), μm1m2=μ(δm1,0δm2,1+δm1,1δm2,0, Γm1m2L(R)=δm1,m2Γm1(L(R), and Um1m2=(δm1,0+δm1,1)δm2,2.

Performing transformation to the field rotational frame [Eq. (18)] with

Ŝ(t)=iω02tn^0n^1
(C1)

leads to a Hamiltonian similar to (19) with ε¯0=ε0+ω0/2, ε¯1=ε1ω0/2, and ε¯2=ε2; with V¯k0(t)=Vk0eiω0t/2, V¯k1(t)=Vk1eiω0t/2, and V¯k2(t)=Vk2; with μE0 inducing electron transfer between orbitals 0 and 1,

12μE0d^1d^0eiϕ0+H.c.,
(C2)

and with coupling to the radiation field—the last line in Eq. (19)—becoming time-dependent,

αVαptd^2d^0eiω0t/2+d^1eiω0t/2âα+H.c..
(C3)

Diagonalizing the 0–1 block of the molecular Hamiltonian yields eigen-orbitals ε+ and ε and leads to expressions for system-baths (contacts and radiation field) couplings in the form

kiVk1sinθei(ω0t+ϕ0)/2iVk0cosθei(ω0t+ϕ0)/2ĉkd^++iVk1cosθei(ω0t+ϕ0)/2+iVk0sinθei(ω0t+ϕ0)/2ĉkd^+Vk2ĉkd^2+H.c.,
(C4)
αiVαptd^2(sinθei(ω0t+ϕ0)/2cosθei(ω0t+ϕ0)/2)d^++  (cosθei(ω0t+ϕ0)/2+sinθei(ω0t+ϕ0)/2)d^+H.c..
(C5)

Analyzing expression (C4) within the RWA-type consideration used in Ref. 64 and neglecting off-diagonal elements in the electronic dissipation rate matrix Γ leads to the diagonal structure of the Green’s function Gm1m2(0)=δm1m2Gm1m1(0) at zero order in the molecule–radiation field coupling.

Thus, the dc part of the absorption spectrum (i.e., dc part of the incoming photon flux) takes the form

2RetdtF<(t,t)G22(0)>(t,t)  ×G++(0)<(t,t)sin2θeiω0(tt)/2+cos2θeıω0(tt)/2  +G(0)<(t,t)cos2θeiω0(tt)/2+sin2θeıω0(tt)/2.
(C6)

Taking integral in t′ leads to Eq. (20).

1.
S. W.
Wu
,
G. V.
Nazin
, and
W.
Ho
,
Phys. Rev. B
77
,
205430
(
2008
).
2.
X. H.
Qiu
,
G. V.
Nazin
, and
W.
Ho
,
Science
299
,
542
(
2003
).
3.
Z.-C.
Dong
,
X.-L.
Guo
,
A. S.
Trifonov
,
P. S.
Dorozhkin
,
K.
Miki
,
K.
Kimura
,
S.
Yokoyama
, and
S.
Mashiko
,
Phys. Rev. Lett.
92
,
086801
(
2004
).
4.
K.
Kimura
,
K.
Miwa
,
H.
Imada
,
M.
Imai-Imada
,
S.
Kawahara
,
J.
Takeya
,
M.
Kawai
,
M.
Galperin
, and
Y.
Kim
,
Nature
570
,
210
(
2019
).
5.
C.
Chen
,
P.
Chu
,
C. A.
Bobisch
,
D. L.
Mills
, and
W.
Ho
,
Phys. Rev. Lett.
105
,
217402
(
2010
).
6.
Y.
Zhang
,
Y.
Luo
,
Y.
Zhang
,
Y.-J.
Yu
,
Y.-M.
Kuang
,
L.
Zhang
,
Q.-S.
Meng
,
Y.
Luo
,
J.-L.
Yang
,
Z.-C.
Dong
, and
J. G.
Hou
,
Nature
531
,
623
(
2016
).
7.
H.
Imada
,
K.
Miwa
,
M.
Imai-Imada
,
S.
Kawahara
,
K.
Kimura
, and
Y.
Kim
,
Nature
538
,
364
(
2016
).
8.
N. L.
Schneider
,
J. T.
,
M.
Brandbyge
, and
R.
Berndt
,
Phys. Rev. Lett.
109
,
186601
(
2012
).
9.
Z.
Ioffe
,
T.
Shamai
,
A.
Ophir
,
G.
Noy
,
I.
Yutsis
,
K.
Kfir
,
O.
Cheshnovsky
, and
Y.
Selzer
,
Nat. Nanotechnol.
3
,
727
(
2008
).
10.
D. R.
Ward
,
N. J.
Halas
,
J. W.
Ciszek
,
J. M.
Tour
,
Y.
Wu
,
P.
Nordlander
, and
D.
Natelson
,
Nano Lett.
8
,
919
(
2008
).
11.
D. R.
Ward
,
D. A.
Corley
,
J. M.
Tour
, and
D.
Natelson
,
Nat. Nanotechnol.
6
,
33
(
2011
).
12.
D.
Natelson
,
Y.
Li
, and
J. B.
Herzog
,
Phys. Chem. Chem. Phys.
15
,
5262
(
2013
).
13.
R.
Chikkaraddy
,
B.
de Nijs
,
F.
Benz
,
S. J.
Barrow
,
O. A.
Scherman
,
E.
Rosta
,
A.
Demetriadou
,
P.
Fox
,
O.
Hess
, and
J. J.
Baumberg
,
Nature
535
,
127
(
2016
).
14.
N.
Kongsuwan
,
A.
Demetriadou
,
R.
Chikkaraddy
,
F.
Benz
,
V. A.
Turek
,
U. F.
Keyser
,
J. J.
Baumberg
, and
O.
Hess
,
ACS Photonics
5
,
186
(
2018
).
15.
M.
Galperin
and
A.
Nitzan
,
Phys. Chem. Chem. Phys.
14
,
9421
(
2012
).
16.
M.
Galperin
,
Chem. Soc. Rev.
46
,
4000
(
2017
).
17.
G. A.
Steele
,
A. K.
Hüttel
,
B.
Witkamp
,
M.
Poot
,
H. B.
Meerwaldt
,
L. P.
Kouwenhoven
, and
H. S. J. v. d.
Zant
,
Science
325
,
1103
(
2009
).
18.
T.
Wagner
,
P.
Talkner
,
J. C.
Bayer
,
E. P.
Rugeramigabo
,
P.
Hänggi
, and
R. J.
Haug
,
Nat. Phys.
15
,
330
(
2019
).
19.
Y.
Selzer
and
U.
Peskin
,
J. Phys. Chem. C
117
,
22369
(
2013
).
20.
M. A.
Ochoa
,
Y.
Selzer
,
U.
Peskin
, and
M.
Galperin
,
J. Phys. Chem. Lett.
6
,
470
(
2015
).
21.
F.
Zhan
,
N.
Li
,
S.
Kohler
, and
P.
Hänggi
,
Phys. Rev. E
80
,
061115
(
2009
).
22.
A.
Donarini
,
A.
Yar
, and
M.
Grifoni
,
Eur. Phys. J. B
85
,
316
(
2012
).
23.
U.
Peskin
,
Fortschr. Phys.
65
,
1600048
(
2017
).
24.
H.
Al Husseini
,
S.
Abdalah
,
K.
Al Naimee
,
R.
Meucci
, and
F.
Arecchi
,
Nanomater. Nanotechnol.
8
,
184798041878238
(
2018
).
25.
V.
Leyton
,
S.
Weiss
, and
M.
Thorwart
,
Phys. Rev. B
97
,
214303
(
2018
).
26.
T.-H.
Park
and
M.
Galperin
,
Europhys. Lett.
95
,
27001
(
2011
).
27.
T.-H.
Park
and
M.
Galperin
,
Phys. Rev. B
84
,
075447
(
2011
).
28.
M. I.
Sena-Junior
,
L. R. F.
Lima
, and
C. H.
Lewenkopf
,
J. Phys. A: Math. Theor.
50
,
435202
(
2017
).
29.
S.
Kohler
,
J.
Lehmann
, and
P.
Hänggi
,
Phys. Rep.
406
,
379
(
2005
).
30.
D.
Thuberg
,
E.
Muñoz
,
S.
Eggert
, and
S. A.
Reyes
,
Phys. Rev. Lett.
119
,
267701
(
2017
).
31.
B. H.
Wu
and
J. C.
Cao
,
Phys. Rev. B
73
,
245412
(
2006
).
32.
T.-S.
Ho
,
K.
Wang
, and
S.-I.
Chu
,
Phys. Rev. A
33
,
1798
(
1986
).
33.
J.
Lehmann
,
S.
Kohler
,
P.
Hänggi
, and
A.
Nitzan
,
Phys. Rev. Lett.
88
,
228305
(
2002
).
34.
J.
Lehmann
,
S.
Kohler
,
P.
Hänggi
, and
A.
Nitzan
,
J. Chem. Phys.
118
,
3283
(
2003
).
35.
B. H.
Wu
and
C.
Timm
,
Phys. Rev. B
81
,
075309
(
2010
).
36.
S.
Kohler
,
J.
Lehmann
, and
P.
Hänggi
,
Superlattices Microstruct.
34
,
419
(
2003
).
37.
S.
Kohler
,
S.
Camalet
,
M.
Strass
,
J.
Lehmann
,
G.-L.
Ingold
, and
P.
Hänggi
,
Chem. Phys.
296
,
243
(
2004
).
38.
G.
Stefanucci
,
S.
Kurth
,
A.
Rubio
, and
E. K. U.
Gross
,
Phys. Rev. B
77
,
075339
(
2008
).
39.
B. H.
Wu
and
J. C.
Cao
,
J. Phys.: Condens. Matter
20
,
085224
(
2008
).
40.
D.
Rai
and
M.
Galperin
,
J. Phys. Chem. C
117
,
13730
(
2013
).
41.
B. H.
Wu
and
J. C.
Cao
,
Phys. Rev. B
81
,
085327
(
2010
).
42.
N.
Tsuji
,
T.
Oka
, and
H.
Aoki
,
Phys. Rev. B
78
,
235124
(
2008
).
43.
A. K.
Eissing
,
V.
Meden
, and
D. M.
Kennes
,
Phys. Rev. B
94
,
245116
(
2016
).
44.
M.
Bukov
,
M.
Kolodrubetz
, and
A.
Polkovnikov
,
Phys. Rev. Lett.
116
,
125301
(
2016
).
45.
I.
Kinoshita
,
A.
Misu
, and
T.
Munakata
,
J. Chem. Phys.
102
,
2970
(
1995
).
46.
M.
Galperin
and
A.
Nitzan
,
J. Chem. Phys.
124
,
234709
(
2006
).
47.
B.
Gu
and
I.
Franco
,
Phys. Rev. A
98
,
063412
(
2018
).
48.
I.
Martin
and
D.
Mozyrsky
,
Phys. Rev. B
71
,
165115
(
2005
).
49.
N.
Sergueev
,
D.
Roubtsov
, and
H.
Guo
,
Phys. Rev. Lett.
95
,
146803
(
2005
).
50.
T.
Frederiksen
,
M.
Paulsson
,
M.
Brandbyge
, and
A.-P.
Jauho
,
Phys. Rev. B
75
,
205413
(
2007
).
51.
R.
Avriller
and
T.
Frederiksen
,
Phys. Rev. B
86
,
155411
(
2012
).
52.
A.
Baratz
,
M.
Galperin
, and
R.
Baer
,
J. Phys. Chem. C
117
,
10257
(
2013
).
53.
H.
Haug
and
A.-P.
Jauho
,
Quantum Kinetics in Transport and Optics of Semiconductors
, second, substantially revised edition (
Springer
,
Berlin Heidelberg
,
2008
), Vol. 123.
54.
G.
Stefanucci
and
R.
van Leeuwen
,
Nonequilibrium Many-Body Theory of Quantum Systems. A Modern Introduction
(
Cambridge University Press
,
2013
).
55.
Y.
Gao
and
M.
Galperin
,
J. Chem. Phys.
144
,
174113
(
2016
).
56.
A.-P.
Jauho
,
N. S.
Wingreen
, and
Y.
Meir
,
Phys. Rev. B
50
,
5528
(
1994
).
57.
P.
Myöhänen
,
A.
Stan
,
G.
Stefanucci
, and
R.
van Leeuwen
,
Phys. Rev. B
80
,
115107
(
2009
).
58.
G.
Baym
and
L. P.
Kadanoff
,
Phys. Rev.
124
,
287
(
1961
).
59.
G.
Baym
,
Phys. Rev.
127
,
1391
(
1962
).
60.
A.
Nitzan
and
M.
Galperin
,
J. Phys. Chem. Lett.
9
,
4886
(
2018
).
61.
S.
Mukamel
and
M.
Galperin
,
J. Phys. Chem. C
123
,
29015
(
2019
).
62.
M.
Galperin
,
M. A.
Ratner
, and
A.
Nitzan
,
J. Chem. Phys.
121
,
11965
(
2004
).
63.
M. O.
Scully
,
Phys. Rev. Lett.
104
,
207701
(
2010
).
64.
M. O.
Scully
,
S.-Y.
Zhu
, and
A.
Gavrielides
,
Phys. Rev. Lett.
62
,
2813
(
1989
).