Optical properties of periodically driven open nonequilibrium quantum systems

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 fluorescence^1–3 and phosphorescence⁴ were reported in the literature. Electroluminescence allowed to visualize intra-molecular interactions,^5,6 explore real space energy transfer,⁷ and probe charge fluctuations in biased molecular junctions.⁸ 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.¹⁷ 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.¹⁸ 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 equation^21–25 and nonequilibrium Green’s function^26–28 approaches. Periodicity in external driving allows mapping of the original time-dependent problem into effective time-independent formulation in an extended Floquet space.²⁹ 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,³⁰ scattering matrix approach,³¹ quantum master equation,^32–35 and Green’s functions^36–40 are available in the literature. For strongly correlated systems, combinations of the Floquet with the slave boson,⁴¹ dynamical mean field theory,⁴² and functional renormalization group⁴³ 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 eV⁴⁵ is of the same order of magnitude or stronger than thermal energy k_BT ∼ 0.01 eV. Thus, the Redfield/Lindblad quantum master equation, which relies on assumption Γ ≪ k_BT, 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 junctions⁴⁶ by accounting for periodic external driving and extends the recent study on optical absorption properties of driven matter⁴⁷ 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.⁴⁸ 

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 usual^49–51 although not always a reliable⁵² level of modeling in first principles simulations. In addition, the molecule is subjected to external driving by monochromatic classical field $E→(t)$⁠,

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.

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

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

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. H^M 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. V_km 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 m₂ to orbital m₁ 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

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

where T_c 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

Here, Hamiltonian H^M, 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,

where

$gk(τ1,τ2)≡−i⟨Tcĉk(τ1)ĉk†(τ2)⟩0$ is the Green’s function of the free electron in state k, and $fα(τ1,τ2)≡−i⟨Tcâα(τ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 I_pt(t)⁵⁵ and electron current I_L(R)(t) through interface L (R),⁵⁶ 

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,

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,⁵⁷ periodic driving allows us to formulate an effective time-independent problem in an extended Floquet space.²⁹ 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

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

Here, $Ga(f;E)≡[Gr(f;E)]†$⁠,

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

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),

where

are the radiative energy dissipation rate and electronic decay rate, respectively, f_L(R)(E) is the Fermi–Dirac distribution, and N_pt(ω) 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 laws^58,59 and may lead to qualitative failures.⁶⁰ Thus, it is not applicable in optical studies of open quantum systems.⁶¹ 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)≡−i∫dEGmm<(f,f;E)/2π$⁠, at subsequent steps, s and s + 1, of the iterative procedure is less than 10⁻⁶: $|nm(s+1)(f)−nm(s)(f)|<10−6$ for each f.

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

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 U₀ and the corresponding unit of flux I₀ = U₀/ℏ. Unless stated otherwise, parameters are the following: k_BT = 0.03, electron escape rates are $ΓmL=ΓmR=0.01$⁠, and the energy dissipation rate is γ = 10⁻⁴. 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. Fermi energy is taken as origin, E_F = 0, and bias V_sd is applied symmetrically μ_L = E_F + |e|V_sd/2 and μ_R = E_F − |e|V_sd/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⁻³.

We start from the model considered in Ref. 47. Here, ε₀ = −0.5, ε₁ = 0.5, and ε₂ = 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.6 [see Fig. 2(a)]. As discussed in Ref. 47, at zero applied bias voltage (V_sd = 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 $μ→⋅E→0∼∑m=12Γm$ (here $Γm≡ΓmL+ΓmR$⁠) and increases with strength of the driving [see Fig. 2(c)].

It is easy to understand the effect by using the rotating frame of the driving field. Because driving is at resonance, ω₀ = ε₂ − ε₁, 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)$⁠,

where $n^m≡d^m†d^m$⁠. The transformation leads to

Here, $ε¯0,1=ε0,1+ω0/2$⁠, $ε¯2=ε2−ω0/2$⁠, $V¯k0(1)(t)=Vk0(1)e−iω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 ω₀/2 for levels 0 and 1 and down by ω₀/2 for level 2. This time-independent problem can be treated within the standard NEGF–SCBA.⁶² 

In the rotating frame of the field, the splitting, $μ→0⋅E→0$⁠, 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 V_sd = 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 $μ→⋅E→0=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.1, ε₁ = 0.1, and ε₂ = 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.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).

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 ε₂ − ε₀ in the absence of driving when the driving is present can have four peaks at frequencies (see  Appendix C for details)

In particular, for the parameters of calculation, one expects three peaks for $μ→⋅E→0=ω0$ and four peaks otherwise [two additional peaks for each $μ→⋅E→0$ 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.5, LUMO ε₁ = 0.4, and LUMO+1 ε₂ = 0.6 levels. External driving couples LUMO and LUMO+1 and is chosen to be in resonance with the optical transition, ω₀ = 0.2. Its coupling strength is taken as $μ→⋅E→0=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 ϕ₀ 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.

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

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

γ(ω) 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

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

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

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

where

$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 $e−iω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 t₁ and t₂, Eq. (B1) can be rewritten as

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

where

Note that in our consideration, t₀ 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 N_pt(ω_p) = 1, Green’s function $Fp<$ for the mode can be written in time domain as

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

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

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)=Vk1e−iω0t/2$⁠, and $V¯k2(t)=Vk2$⁠; with $μ→⋅E→0$ inducing electron transfer between orbitals 0 and 1,

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

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