CHEERS: A tool for Correlated Hole-Electron Evolution from Real-time Simulations
CCHEERS: A tool for Correlated Hole-Electron Evolution from Real-time Simulations
E. Perfetto
1, 2 and G. Stefanucci
2, 3 CNR-ISM, Division of Ultrafast Processes in Materials (FLASHit),Area della ricerca di Roma 1, Monterotondo Scalo, Italy Dipartimento di Fisica, Universit`a di Roma Tor Vergata,Via della Ricerca Scientifica, 00133 Rome, Italy INFN, Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma, Italy (Dated: October 9, 2018)We put forward a practical nonequilibrium Green’s function (NEGF) scheme to perform real-timeevolutions of many-body interacting systems driven out of equilibrium by external fields. CHEERSis a computational tool to solve the NEGF equation of motion in the so called generalized Kadanoff-Baym ansatz and it can be used for model systems as well as first-principles Hamiltonians. Dynam-ical correlation (or memory) effects are added to the Hartree-Fock dynamics through a many-bodyself-energy. Applications to time-dependent quantum transport, time-resolved photoabsorption andother ultrafast phenomena are discussed.
I. INTRODUCTION
Although the laws of quantum mechanics have beenformulated almost a century ago, the behavior of quan-tum matter under (extreme) nonequilibrium conditionsremains still largely unexplored. Modern advances inlaser technology [1–3] make today possible to film, withan unprecedented time resolution, the genesis and devel-opment of photoemission processes, exciton formation,charge transfers, charge migrations, Auger decays andother ultrafast phenomena. This is the realm of attosec-ond physics which calls for accurate theories and efficientnumerical schemes to predict the evolution of many-bodyquantum systems.One of the most versatile formalism to deal withthe quantum many-body problem is the diagrammaticGreen’s function theory. The extension to out-of-equilibrium situations is known as the Non-EquilibriumGreen’s Function (NEGF) theory [4–6] and the funda-mental equations, known as Kadanoff-Baym equations(KBE), date back to the mid sixties [7–9]. Despite theenormous advance in computational capabilities the KBEare still rather burdensome to solve numerically. Infact, their implementations have been so far restrictedto atoms, diatomic molecules or model systems [10–19].In the mid-1980s Lipavsky et al. [20] proposed the socalled Generalized Kadanoff-Baym Ansatz (GKBA) tocollapse the KBE for the two-times Green’s function intoa single equation for the one-time one-particle densitymatrix, thus drastically reducing the computational cost.The GKBA is exact in the Hartree-Fock (HF) approxima-tion and it is expected to be accurate when the averagetime between two consecutive collisions is longer thanthe quasiparticle decay time (see Ref. [21] for a recentdiscussion).The appealing feature of the GKBA is that the NEGFformalism is converted into a time-dependent density-matrix functional theory [22–28] which shares a funda-mental property with many-body perturbation theory,i.e., the systematic inclusion of correlations through a proper selection of self-energy diagrams. Recent ap-plications of the NEGF+GKBA approach include thenonequilibrium dynamics [29, 30] and many-body local-ization [31] of Hubbard clusters, equilibrium absorptionof sodium clusters [32], transient absorption [33–35] andcarrier dynamics [36, 37] of semiconductors.In this work we describe CHEERS, a first-principlesnumerical tool based on NEGF+GKBA to simulatethe time evolution of interacting systems. CHEERStime-evolutions contain dynamical correlation (or mem-ory) effects responsible for double (and multiple) excita-tions, decoherence-induced charge separation, Auger de-cays, shake-up dynamics, image-charge renormalizations,etc. Standard time-dependent HF simulations are recov-ered by simply switching off the effects of correlations.CHEERS has been already used to study the charge dy-namics of model molecular junctions [21], the transientphotoabsorption spectrum of noble gas atoms [38], theformation of charge-transfer excitons and their subse-quent separation in donor-acceptor complexes [39], theattosecond pulse-induced charge migration in the pheny-lalanine aminoacid [40] and time-resolved Auger de-cays [41].The paper is organized as follows. In Section II we dis-cuss the physical systems of interest and write down themany-body Hamiltonian to describe them. The theoret-ical framework based on NEGF and GKBA is outlinedin Section III along with the NEGF equation of motionsolved by CHEERS. In Section IV we list the observ-able quantities accessible from the solution of the NEGFequation. A description of the implementation details isgiven in Section V. Conclusions and outlooks are drawnin Section VI.
II. PHYSICAL SYSTEMS
We consider a quantum system with a finite numberof nuclei. For the time being we fix the nuclear coordi-nates and focus on the electronic degrees of freedom. Let { ϕ i ( x ) } be a set of localized orthonormal spin-orbitals a r X i v : . [ c ond - m a t . o t h e r] O c t suited to describe the equilibrium and nonequilibriumproperties of the bound electrons. In our notation x = r σ comprises a spatial coordinate r and a spin-projection σ ,and ϕ i ( x ) = (cid:104) i | x (cid:105) . The localized states could be gen-erated by orthonormalizing a set of Slater-type orbitals(STO) or Gaussian-type orbitals (GTO) centered aroundeach nucleus or they could be the bound Hartree-Fock(HF) orbitals or the bound Kohn-Sham (KS) orbitals re-sulting from some self-consistent HF or KS calculationrespectively. We will give more details on the possiblechoice of the basis set in Section V. For each spin-orbital ϕ i we define the corresponding annihilation and creationoperator ˆ c i and ˆ c † i . The equilibrium Hamiltonian of thesystem in second quantization then readsˆ H eqsys = (cid:88) ij h eq ij ˆ c † i ˆ c j + 12 (cid:88) ijmn v ijmn ˆ c † i ˆ c † j ˆ c m ˆ c n . (1)Here h eq ij are the matrix elements of the single-particleHamiltonian (atomic units are used throughout): h eq ij ≡ (cid:104) i | ˆ p V n + ˆ V SO | j (cid:105) , (2)with ˆ V n the nuclear potential and ˆ V SO the spin-orbit in-teraction potential. The second term in Eq. (1) describesthe electron-electron interaction with Coulomb integrals v ijmn ≡ (cid:90) d x d x (cid:48) ϕ ∗ i ( x ) ϕ ∗ j ( x (cid:48) ) ϕ m ( x (cid:48) ) ϕ n ( x ) | r − r (cid:48) | , (3)where (cid:82) d x = (cid:82) d r (cid:80) σ .We are interested in studying the quantum evolutioninduced by an external electromagnetic field with spatialvariations on length-scales much longer than the lineardimension of the system. For nanometer-sized moleculesthis condition implies photon energies up to 1 keV andhence the possibility of photoionization. To describe pho-toionization processes it is necessary to extend the local-ized basis and include delocalized spin-orbitals { ϕ µ ( x ) } for electrons in the continuum. Without loss of general-ity, we choose the ϕ µ ’s as eigenstates of the free-particleHamiltonian far away from the system boundaries andwe denote by (cid:15) µ their energy. We also require that the ϕ µ ’s are orthogonal to the ϕ i ’s and orthonormal betweenthemeselves, i.e., (cid:104) µ | µ (cid:48) (cid:105) = δ µµ (cid:48) . An example of how toconstruct the ϕ µ ’s is given in Section V. We discard theCoulomb interaction between two electrons in the contin-uum since for weak pulses double ionization is stronglysuppressed. Of course, if X-rays are used then a sec-ond (Auger) electron can be ejected. As we discuss inSection III, in this case the approximation is justifiedprovided that the photoelectron and the Auger electronhave different energies. The continuum Hamiltonian thenreads ˆ H cont = (cid:88) µ (cid:15) µ ˆ c † µ ˆ c µ . (4) Taking into account that the electric field E is uniform forall localized states, we can write the interaction Hamil-tonian between light and matter asˆ H E ( t ) = ˆ H E sys ( t ) + ˆ H E ion ( t ) , (5)where ˆ H E sys ( t ) = E ( t ) · (cid:88) ij d ij ˆ c † i ˆ c j (6)is the part responsible for reshuffling the electrons be-tween localized states whereasˆ H E ion ( t ) = E ( t ) · (cid:88) iµ (cid:16) d iµ ˆ c † i ˆ c µ + h . c . (cid:17) (7)is the part responsible for photoionization processes. InEqs. (6) and (7) the vector of matrices d is the dipolemoment defined as d ab = (cid:90) d x ϕ ∗ a ( x ) r ϕ b ( x ) , (8)where a and b are indices either in the localized { i } ordelocalized { µ } sector. In Eq. (5) we are discarding tran-sitions µ → µ (cid:48) between delocalized electrons since we aremainly concerned with ultrafast fields (by the time thepopulation of a µ -state becomes relevant the electric fieldvanishes).During the first few femtoseconds after ionization froma (semi) core state, the Auger decay is one of the most rel-evant recombination channels. To account for the Augereffect in our description (and hence to deal with soft X-ray pulses too) we include the Coulomb matrix elementsresponsible for two localized (valence) electrons to scat-ter in one localized (core) electron and one delocalized(continuum) electron. In second quantization the Augerinteraction Hamiltonian readsˆ H Auger = (cid:88) ijm (cid:88) µ v Aijmµ (cid:16) ˆ c † i ˆ c † j ˆ c m ˆ c µ + h . c . (cid:17) , (9)where the Coulomb integrals v Aijmµ are defined as inEq. (3) with ϕ n → ϕ µ .In CHEERS the quantum system can also be contactedto metallic leads with which to exchange electrons and en-ergy. Arbitrary time-dependent voltages can be appliedto the leads to study transient currents, steady-states,AC responses or other kind of transport properties. It isalso possible to switch on a thermomechanical field [42] tocalculate time-dependent thermal currents [43–45]. Thecombination of applied voltages and external laser pulses E ( t ) can instead be used to access the optical propertiesof current-carrying molecular junctions [46]. In all casesthe leads are treated as noninteracting semi-infinite crys-tals with a finite cross section and described in termsof semi-infinite Bloch states ϕ αk ( x ), where α is the leadindex and k specifies the energy (cid:15) αk of the Bloch state.The Bloch states have to be orthogonal to both the ϕ i ’s i j mn µ Augerelectron ⌫ ↵k E ( t ) photoelectronconductingelectron V ↵ ( t ) T i,↵k v Aijmµ v ijmn FIG. 1: Illustration of the possible physical effects that canbe addressed with the Hamiltonian in Eq. (12). and the ϕ µ ’s. Thus, in a quantum transport setup the µ -states are the free-particle continuum states of the quan-tum system contacted to leads. The second-quantizedform of the leads Hamiltonian isˆ H lead ( t ) = (cid:88) α (cid:88) k Φ α ( t ) ( (cid:15) αk + V α ( t )) ˆ c † αk ˆ c αk , (10)where V α ( t ) is the applied voltage and Φ α ( t ) = T α ( t ) /T is the ratio between the temperature at time t and theequilibrium temperature (the thermomechanical field istherefore Ψ α = Φ α − H tun = (cid:88) i,αk (cid:16) T i,αk ˆ c † i ˆ c αk + h . c . (cid:17) , (11)where T i,αk is the tunneling amplitude, i.e., the matrixelement of the single-particle Hamiltonian, between thestates ϕ i and ϕ αk .To summarize, the full Hamiltonian has the formˆ H = ˆ H sys + ˆ H cont + ˆ H lead + ˆ H E ion + ˆ H Auger + ˆ H tun , (12)where ˆ H sys = ˆ H eqsys + ˆ H E sys . (13)In Fig. 1 we show an illustration of the physics that canbe studied with the Hamiltonian in Eq. (12). III. GREEN’S FUNCTION FORMULATION
Although the full Hamiltonian in Eq. (12) includes onlythe Auger scattering between bound electrons and con-tinuum electrons, neglects the electron-electron interac-tion in the leads and discards the effects of E ( t ) on thefree-particle Hamiltonian of the continuum states, thegeneral solution of the problem is still challenging. Tomake some progress we restrict the physical situations ofinterest. In experiments the momentum of the photoelec-tron is often well separated from the momentum of theAuger electron. Having in mind this type of experiments we take the set S ion of photoelectron states and the set S Auger of Auger-electron states as two disjoint sets. Thisimplies that in Eq. (7) we can restrict the sum over µ to states with µ ∈ S ion and, similarly, in Eq. (9) we canrestrict the sum over µ to states with µ ∈ S Auger . A. NEGF equations
The method of choice to investigate the electron dy-namics is the Non-Equilibrium Green’s Function (NEGF)approach [4–6]. The Green’s function G ij ( z, z (cid:48) ) withtimes z, z (cid:48) on the Keldysh contour and indices i, j in thelocalized sector satisfies the equation of motion (cid:20) i ddz − h HF ( z ) (cid:21) G ( z, z (cid:48) ) = δ ( z, z (cid:48) ) + I emb ( z, z (cid:48) )+ I ion ( z, z (cid:48) ) + I coll ( z, z (cid:48) ) + I Auger ( z, z (cid:48) ) (14)and its adjoint. Let us discuss the quantities in Eq. (14).The matrix h HF is the single-particle Hartree-Fock (HF)Hamiltonian with elements h HF ,ij ( z ) ≡ h eq ij + V HF ,ij ( z ) + E ( z ) · d ij , (15)where the HF potential V HF ,ij ( z ) ≡ (cid:88) mn ρ nm ( z ) w imnj ( z ) (16)is expressed in terms of the one-particle density matrix ρ nm ( z ) ≡ − iG nm ( z, z + ) (17)and the difference between the direct and exchangeCoulomb integrals w imnj ( z ) ≡ v imnj ( z ) − v imjn ( z ) . (18)Hereafter we consider the more general case of a time-dependent interaction v = v ( z ), useful to deal with adia-batic switching, interaction quenches, etc..The embedding / ionization integral I emb / ion ( z, z (cid:48) ) ≡ (cid:90) d ¯ z Σ emb / ion ( z, ¯ z ) G (¯ z, z (cid:48) ) (19)is a convolution on the Keldsyh contour between the em-bedding/ionization self-energy and the Green’s function.The embedding self-energy is responsible for tunneling ofelectrons to/from the leads and reads [14, 47, 48]Σ emb ,ij ( z, ¯ z ) = (cid:88) αk T i,αk ( z ) g αk ( z, ¯ z ) T αk,j (¯ z ) , (20)where we allow for a time-dependent tunneling amplitude T i,αk ( z ). The ionization self-energy is instead responsiblefor the photoionization of the system and reads [38, 40]Σ ion ,ij ( z, ¯ z ) = (cid:88) µ ∈ S ion ( E ( z ) · d iµ ) g µ ( z, ¯ z ) ( E (¯ z ) · d µj ) . (21) = + ij i n m jp q r s G nm G pq G sr v i r pn v m q s j i n m jq s r p v irpn v mqsj G pm G nq G sr Σ FIG. 2: Diagrammatic representation of the 2B self-energy.Wiggly lines denote the Coulomb interaction v . Both self-energies are expressed in terms of a free-particleGreen’s function g . For the embedding self-energy g αk isthe solution of the equation of motion (cid:20) i ddz − Φ α ( z ) ( (cid:15) αk − V α ( z )) (cid:21) g αk ( z, ¯ z ) = δ ( z, ¯ z ) , (22)whereas for the ionization self-energy g µ is the solutionof the equation of motion (cid:20) i ddz − (cid:15) µ (cid:21) g µ ( z, ¯ z ) = δ ( z, ¯ z ) . (23)All equations of motion, including Eq. (14), are solvedwith the appropriate Kubo-Martin-Schwinger boundaryconditions [5].There are two more terms to be discussed. The firstterm is the collision integral I coll ( z, z (cid:48) ) ≡ (cid:90) d ¯ z Σ( z, ¯ z ) G (¯ z, z (cid:48) ) , (24)where the correlation self-energy Σ = Σ[ v, G ] is a func-tional of the interaction v and the Green’s function G .The exact Σ is the sum of all skeletonic self-energy dia-grams with propagators G and interaction lines v [5]. InCHEERS this self-energy is implemented at the level ofthe second-Born (2B) approximation, i.e.,Σ ij ( z, ¯ z ) = (cid:88) mn pq sr G mn ( z, ¯ z ) G pq ( z, ¯ z ) G sr (¯ z, z ) × v irpm ( z ) w nqsj (¯ z ) . (25)In Fig. 2 we show the corresponding diagrammatic rep-resentation. We mention that the 2B approximation hasbeen successfully applied to equilibrium spectral prop-erties [49] and total energies [50] of molecular systems.Comparisons against numerically accurate real-time sim-ulations of 1D systems [18, 51] and weakly correlatedmodel nanostructures of different geometries [29, 30, 39,52–56] indicate that the 2B approximation remains accu-rate even out of equilibrium.The last term I Auger can be written as in Eq. (24)but the self-energy contains all diagrams with at leastone Auger interaction line v A [57, 58]. If we are onlyinterested in describing the Auger physics then we canapproximate I Auger ( z, ¯ z ) = (cid:90) d ¯ z Σ Auger ( z, ¯ z ) G (¯ z, z (cid:48) ) , (26) where the Auger self-energy reads [41]Σ Auger ,ij ( z, ¯ z ) = (cid:88) mn pq (cid:88) µν ∈ S Auger G mn ( z, ¯ z ) × (cid:2) G µν ( z, ¯ z ) G pq (¯ z, z )( v Aiqmµ w Aνnpj + v Aiqµm w Anνpj )+ G pq ( z, ¯ z ) G µν (¯ z, z ) v Aiνpm w Anqµj (cid:3) . (27)This self-energy corresponds to the 2B approximationwith interaction lines v A and follows from the many-body identity Σ = − iv A G G − where the two-particleGreen’s function G describes a single Auger scattering.Notice that for strong valence-valence repulsive energies(typically well above 1 eV) it is crucial to replace the(first-order) single-scattering approximation to G withthe T -matrix approximation in the particle-particle chan-nel [59, 60]. Currently, CHEERS does not contain im-plementations of the T -matrix approximation and henceit can only be used to study Auger decays in molecules(e.g., organic molecules) with fairly delocalized valenceorbitals .We observe that Σ Auger depends on one Green’s func-tion with both indices in the continuum. ThereforeEq. (14) for G ij can be solved only provided that wecouple it to an equation for G µν . To second order in v A the equation of motion for G µν with µ, ν ∈ S Auger is (cid:20) i ddz − (cid:15) µ (cid:21) G µν ( z, z (cid:48) ) = δ ( z, z (cid:48) ) + (cid:88) ρ ∈ S Auger (cid:88) mn pq sr × (cid:90) d ¯ z G mn ( z, ¯ z ) G pq ( z, ¯ z ) G sr (¯ z, z ) v Aµrpm w Anqsρ × G ρν (¯ z, z ) (28)and it is linear in G µν . An important consequence of thisresult is that if the Green’s function of the initial state isblock-diagonal in the indices i and µ , i.e., G iµ = 0, thenit remains block diagonal.Equations (14) and (28) form a coupled system of non-linear integro-differential equations and except for the 2Bself-energy approximation no other approximations havebeen made. The full numerical solution of these equa-tions requires to convert the contour-time G to real-time G ’s, a procedure leading to the so-called Kadanoff-Baymequations (KBE). The numerical solution of the KBEis rather demanding, especially for large basis sets andmany nonvanishing four-index Coulomb integrals. In thenext section we discuss how the computational cost isdrastically reduced by the GKBA. B. GKBA equations
We rewrite the equation of motion (14) as (cid:20) i ddz − h HF ( z ) (cid:21) G ( z, z (cid:48) ) = δ ( z, z (cid:48) )+ (cid:90) d ¯ z Σ tot ( z, ¯ z ) G (¯ z, z (cid:48) ) , (29)with total self-energyΣ tot = Σ ion + Σ emb + Σ + Σ Auger . (30)In Σ tot only the last two terms are functionals of G .Choosing z and z (cid:48) on different branches of the Keldyshcontour we obtain the KBE for the lesser and greaterGreen’s functions [4–6] (cid:20) i ddt − h HF ( t ) (cid:21) G ≶ ( t, t (cid:48) ) = (cid:90) d ¯ t Σ Rtot ( t, ¯ t ) G ≶ (¯ t, t (cid:48) )+ (cid:90) d ¯ t Σ ≶ tot ( t, ¯ t ) G A (¯ t, t (cid:48) ) , (31)where, for any function F , the superscript R and Adenotes the retarded and advanced components respec-tively: F R / A ( t, t (cid:48) ) = ± θ ( ± t ∓ t (cid:48) )[ F > ( t, t (cid:48) ) − F < ( t, t (cid:48) )] . (32)Similarly, from the adjoint of Eq. (14) we find G ≶ ( t, t (cid:48) ) (cid:34) i ←− ddt (cid:48) − h HF ( t (cid:48) ) (cid:35) = (cid:90) d ¯ t G R ( t, ¯ t )Σ ≶ tot (¯ t, t (cid:48) )+ (cid:90) d ¯ t G ≶ ( t, ¯ t )Σ Atot (¯ t, t (cid:48) ) . (33)Without any loss of generality we assume that the sys-tem is in equilibrium until a certain time t switch > V α ( t ) = E ( t ) = 0 for t < t switch . To obtain the correlated and contacted (to leads, if any) Green’s func-tion we solve Eqs. (31) and (33) with: • initial condition G ≶ (0 ,
0) given by the HF (henceuncorrelated) lesser/greater Green’s function of theuncontacted system • self-energies calculated using a TD interaction v ( t ) = s ( t ) v and tunneling amplitude T ( t ) = s ( t ) T where s ( t ) is a slow and smooth switching functionbetween the times t = 0 and t = t switch , typically s ( t ) = sin ( πt t switch ).The time t switch is therefore a convergence parameter tobe chosen in such a way that the observables of interestare constant in the absence of external fields for times t > t switch . This initial time-propagation serves to buildup correlations in the inital state. Since v = T = 0 fortimes t < ∞ , i.e., (cid:82) d ¯ t = (cid:82) ∞ d ¯ t .Subtracting Eq. (33) from Eq. (31) and setting t (cid:48) = t we obtain the equation of motion for the one-particledensity matrix ρ ij = − iG
From the solution of Eqs. (43) we can calculate severalquantities of physical interest. The most straighforwardones are the local density n ( r , t ) = (cid:88) ij (cid:88) σ ϕ ∗ i ( r σ ) ϕ j ( r σ ) ρ ji ( t ) , (44)spin density s ( r , t ) = (cid:88) ij (cid:88) σσ (cid:48) ϕ ∗ i ( r σ ) σ σσ (cid:48) ϕ j ( r σ (cid:48) ) ρ ji ( t ) (45)local paramagnetic current j ( r , t ) = 12 (cid:88) ij (cid:88) σ Im [ ϕ ∗ i ( r σ ) ∇ ϕ j ( r σ ) ρ ji ( t )] , (46)and, more generally, any one-body observable. In Fig. 3we show the snapshots of the density variation in thephenylalanine aminoacid induced by an ionizing attosec-ond XUV pulse [40], see Section V B for the implemen-tation details. FIG. 3: Snapshots of the density variation in the phenylala-nine aminoacid induced by an ionizing attosecond XUV pulse.The excess of hole density (blue) and electron density (red)refer to the density averaged over the full time simulation.Reprinted figure with permission from [40]. Copyright 2018by the American Chemical Society.
Depending on the physical problem these basic quan-tities can be further manipulated to calculate typical ex-perimental outcomes. In the following we discuss someof them.
A. Quantum Transport
In molecular electronics the Hamiltonian of the systemdescribes a junction connecting two leads (a source anda drain) kept at a potential difference V . One is usuallyinterested in the total current I flowing through a surface S perpendicular to the electron stream I ( t ) = (cid:90) S d r j (cid:107) ( r , t ) , (47)with j (cid:107) the longitudinal component of j , see Eq. (46).For a DC bias I ( t ) attains a steady value as t → ∞ andthis value can be used to calculate the I − V charac-teristics or the differential conductance G = dI/dV . Ofcourse, the full time-evolution provides other useful infor-mation. The characteristic time to reach a steady stateand the frequencies of the transient oscillations are justtwo examples. Controlling these properties is crucial toengineer ultrafast molecular devices.Time-dependent potential differences V ( t ) do not bringadditional complications nor an increased computationaleffort. We can, for instance, superimpose an AC volt-age of frequency Ω to a DC voltage, i.e., V ( t ) = V DC + V AC sin(Ω t ), and calculate the averaged current as wellas its Fourier coefficients as functions of V DC , V AC andΩ. In addition to provide an alternative to Floquetschemes [61], working in the time domain is particularlyadvantageous to deal with systems perturbed by multi-chromatic drivings. This is the case of, e.g., AC trans-port with superconducting leads [62] or optical spectraof junctions under AC voltages [46].The current I ( t ) can be evaluated either at an inter-face passing through the junction, in accordance withEq. (47), or at the interface with the α lead through theMeir-Wingreen formula [47, 48] I α ( t ) = 4Re (cid:90) t d ¯ t Tr (cid:2) Σ >α ( t, ¯ t ) G < (¯ t, t ) − Σ <α ( t, ¯ t ) G > (¯ t, t ) (cid:3) (48)where Σ α is the α -th contribution to the embedding self-energy of Eq. (20) and G ≶ are calculated from ρ throughthe GKBA. In Ref. 21 we showed that the GKBA resultsfor the current at the interfaces are in excellent agreementwith the full KBE results provided that the bias differ-ence is much smaller than the bandwidth of the leads.We mention that an improved version of the GKBA hasbeen recently proposed to deal with nontrivial spectralstructures of the leads density-of-states [63], thus widen-ing the potential applications of the GKBA in quantumtransport.The energy current J α ( t ), defined as the rate of changeof the energy of lead α , can be calculated similarly to thecharge current in Eq. (49). It is a matter of simple algebrato show that [43] J α ( t ) = 4Im (cid:90) t d ¯ t Tr (cid:104) ˙Σ >α ( t, ¯ t ) G < (¯ t, t ) − ˙Σ <α ( t, ¯ t ) G > (¯ t, t ) (cid:105) (49)where ˙Σ ≶ α ( t, ¯ t ) ≡ ddt Σ ≶ α ( t, ¯ t ). B. Transient Photoabsorption Spectra
In a transient photoabsorption spectrum the systemis driven out of equilibrium by a strong laser pulse (thepump) and successively the intensity per unit frequencyof the transmitted light of a second weak pulse (theprobe) is measured. The resulting transient spectrum de-pends on the shape and duration of the pump and probepulses as well as from the delay τ between these pulses.It is therefore clear that the theoretical calculation of atransient spectrum calls for a time-dependent approach.Let (cid:104) δ d ( t, τ ) (cid:105) be the change of the dipole moment ofthe pump-driven system induced by an electric probefield e ( t ) impinging the system with a delay τ with re-spect to the pump. Then the transient spectrum is givenby S ( ω, τ ) = − (cid:104) ω ˜ e ∗ ( ω ) · (cid:104) δ ˜ d ( ω, τ ) (cid:105) (cid:105) (50)where we have used the convention that quantities witha tilde denote the Fourier transform of the correspondingtime-dependent quantities.The dipole moment of the system can easily becalculated from the one-particle density matrix as ω ( e V ) τ (f s ) HF Kr ω ( e V ) τ (f s ) Kr Kr FIG. 4: Transient photoabsorption spectrum (normalized tothe maximum height) of a krypton gas in the HF (top panel)and 2B (bottom panel) approximation. Reprinted figure withpermission from [38]. Copyright 2015 by the American Phys-ical Society. (cid:104) d ( t ) (cid:105) = (cid:80) ij d ij ρ ji ( t ). The probe-induced dipole mo-ment (cid:104) δ d ( t, τ ) (cid:105) is therefore the difference between thedipole moment generated by a simulation with pump andprobe and the dipole moment generated by a simulationwith only the pump. We observe that the pump pulse caneither bring the system in an excited state of bound elec-trons or generate a multiply ionized system through theionization integral I ion . In the latter case the transientspectrum reveals features about the initial dynamics ofthe expelled photoelectrons [64]. In Fig. 4 we display HF(top) and correlated (bottom) calculations of the tran-sient spectrum of a gas of Kr atoms initially ionized byfew-cycle NIR pump and subsequently probed by an at-tosecond XUV pulse at different delays τ [38], see Sec-tion V A for the implementation details. In HF no signof multiple ionization is visible. On the contrary, thecorrelated 2B results clearly show the absorption line ofKr ions raising up a few femtosecond later than theabsorption line of Kr . C. Transient Photocurrent and Auger Current
In the absence of leads (closed system) the total electriccurrent flowing out of the system, i.e., in the continuumstates, can be calculated from the rate of change of thetotal number of particles in the system I ion ( t ) = dN ( t ) dt = ddt Tr[ ρ ( t )] . (51)If we are interested in resolving the photocurrent accord-ing to the energy of the photoelectrons (photoemissionspectrum) we should include explicitly the continuumstates µ ∈ S ion in the simulation (instead of using Σ ion ).We would then get a third equation for the occupations f µ with µ ∈ S ion coupled to Eqs. (43). Of course thisprocedure is feasible only provided that the energy win-dow of the photoelectrons is not too wide and that thenumber of µ states for the required energy resolution isnot too large. At present the explicit inclusion of photo-electron states has been tested only in one-dimensionalmodel systems [41].The Auger self-energy accounts for processes where twovalence electrons scatter and, after the scattering, endup in a core state and in a continuum state. The rateof growth of the occupation of the continuum state µ ∈S Auger defines the Auger current I Auger ( t ) = ddt f µ ( t ) . (52)Through knowledge of ρ ij and f µ we can follow in realtime the Auger scattering and extract useful informationabout the Auger process, e.g., core relaxation time, shapeof the outgoing density wavepacket, rearrangement of thecore-excited system, etc [41]. V. IMPLEMENTATION DETAILS
From Eqs. (43) and the definition of the various quan-tities therein we find useful to split the input parametersinto five different groups • System : matrix elements h eq ij , d ij and Coulomb in-tegrals v ijmn • Ionization : matrix elements d iµ and energies (cid:15) µ with µ ∈ S ion • Auger : matrix elements d iµ , energies (cid:15) µ andCoulomb integrals v Aijmµ with µ ∈ S Auger • Leads : tunneling amplitudes T i,αk and energies (cid:15) αk • External fields : electric field E ( t ), bias V α ( t ) andtemperature T α ( t ).Currently, real-time simulations in the presence ofleads (open systems) are performed only for modelHamiltonians. Here, the input parameters are set man-ually and can be varied at will. Real-time simulationsbased on first-principles input parameters are possiblefor closed systems (Σ emb = 0) like atoms and moleculesin external laser fields. CHEERS handles the input pa-rameters in different way depending on the nature of thesingle-particle basis set. In the following two subsectionswe describe how CHEERS processes the input generatedin a basis of localized orbitals like, e.g., Slater Type Or-bitals (STO) or Gaussian Type Orbitals (GTO), and ina basis of Kohn-Sham orbitals. In all cases the first stepof CHEERS is to obtain the information contained in System , Ionization , Auger and
Leads . With this informa-tion CHEERS calculates all self-energies and then passesthem to the time-propagation routine, see blue arrows inFig. 5.As illustrated in the left red box of Fig. 5 the time-propagation routine needs also other information. Twomain flags specify the type of evolution. One flag es-tablishes the level of correlation: it can be either aHF evolution – with the extra option of using the HFHamiltonian h HF [ ρ ( t switch )]( t ) with frozen ρ instead of h HF [ ρ ( t )]( t ) – or a correlated 2B evolution. Another flagsets which propagator G R [ ρ ] is used in the GKBA, seesubsection V D for the possible choices. Prior to the timeevolution we also specify a few convergence parameters.The most important ones are the time-step, the switchingtime t switch for the correlation build-up and the number ofpredictor correctors for each time step. Finally, we spec-ify the driving fields in the input External fields . Thereare no restrictions on the time-dependent functions E ( t ), V α ( t ) and T α ( t ), and the computational effort does notdepend on the choice of these functions. During the timestepping the density matrix ρ ij and f µ are either saved orprocessed to generate the output described in Section IV. A. Localized basis
For a description in terms of N single-particle local-ized states {| i (cid:105)} like, e.g., the STO or GTO states,CHEERS needs the matrix elements of the equilibriumHamiltonian h eq ij , dipole vector d ij , overlap matrix S ij = (cid:104) i | j (cid:105) and Coulomb integrals v ijmn . The first step ofCHEERS is to orthonormalize the basis according to | i (cid:105) → (cid:80) m | m (cid:105) S − / mi , calculate h eq ij , d ij and v ijmn in thenew orthonormal basis and run a HF self-consistent cal-culation. The second step is to calculate h eq ij , d ij and v ijmn in the HF basis which, by definition, diagonalizesthe HF Hamiltonian h HF ,ij = (cid:15) HF i δ ij . This is done onlyfor HF states with energy Λ min < (cid:15) HF i < Λ max whereΛ min and Λ max are two convergence cutoff parameters.Of course if Λ min is smaller than the minimum HF eigen- Local basis (STO, GTO) KS basis Models HF States assignment States assignment ⌃ ion ⌃ Auger ⌃ emb (exact or time local) States assignment { c } , { i } , { µ ion } { µ Auger } { c } , { i } , { µ ion } (time local) ⌃ ion ⌃ Auger = 0 ⌃ emb = 0 ⌃ ion ⌃ Auger ⌃ emb = 0 TIME PROPAGATION : frozen HF, HF, 2B: choice of propagator: time-step, switching-time, nr. predictors { µ Auger } { c } , { i } , { µ ion } : laser, bias, temperature OUTPUT
FIG. 5: Architecture of CHEERS. The input contained in
System , Ionization , Auger and
Leads can be either generated fromfirst-principles calculations in a localized basis (left) or KS basis (middle) or, alternatively, it can be set manually for modelsystem calculations (right). The level of correlation, type of propagator, convergence parameters and the input of
Externalfields is given in the red box (bottom right). During the time propagation the density matrix is either saved or processed togenerate the output of interest. value and Λ max is larger than the maximum HF eigen-value then all N HF states are included. On the con-trary, states i = c of energy (cid:15) HF c < Λ min are treated ascore states, i.e., ρ cc (cid:48) = δ cc (cid:48) . Then CHEERS calculatesonly v iccj and v icjc and adds to the equilibrium Hamilto-nian the HF potential generated by the frozen core (fc)electrons: h eq ij → h eq+fc ij = h eq ij + (cid:88) c ( v iccj − v icjc ) . (53)The HF states i = µ with energy (cid:15) HF µ > Λ max are treatedas noninteracting and considered as states of the contin-uum. Accordingly, CHEERS calculates only the dipolematrix elements d iµ later used to construct the ioniza-tion self-energy of Eq. (21). The separation of the HFstates is illustrated in the left green box of Fig. 5.After this preliminary treatment the size of the one-particle density matrix ρ in Eq. (34) becomes N bound = N − N c − N ion , where N c is the number of HF core statesand N ion is the number of µ -states. The STO or GTOdescription of the continuum is, in general, too poor for simulating Auger scattering processes which are thereforenot included, see left blue box of Fig. 5. This meansthat for a localized basis CHEERS solves only the firstof Eqs. (43). For problems involving Auger scattering seenext Section.The µ -states of the localized basis are used to constructan approximate ionization self-energy, thus accountingfor possible photoelectrons due to external laser fields.Of course, if the expected energy range of the photoelec-trons is not covered by the (cid:15) HF µ then the ionization rateis severely underestimated. Meaningful real-time simu-lations do therefore require that the laser frequency isat least smaller than max { (cid:15) HF µ } − π/T where T is theduration of the ionizing pulse.Let us discuss the ionization self-energy. In generalthe STO or GTO basis returns only a few continuumstates; consequently, a photoelectron would be soon re-flected back. This difficulty can be overcome providedthat the ionizing laser pulse is well centered around somefrequency ω P . From Eq. (21) the lesser part of Σ ion van-0ishes whereas the greater part is given byΣ > ion ( t, t (cid:48) ) = (cid:88) ab E a ( t ) σ ab ( t − t (cid:48) ) E b ( t (cid:48) ) , (54)where E a is the a -th component of the electric field E =( E x , E y , E z ) and the tensor σ abij ( t − t (cid:48) ) ≡ − i (cid:88) µ ∈ S ion d aiµ e − i(cid:15) HF µ ( t − t (cid:48) ) d bµj , (55)depends exclusively on the matrix elements of the com-ponents d a of the dipole moment d = ( d x , d y , d z ). TheFourier transform of σ abij reads˜ σ abij ( ω ) = − πi (cid:88) µ ∈ S ion d aiµ δ ( ω − (cid:15) HF µ ) d bµj ≈ i (cid:88) µ ∈ S ion d aiµ Im (cid:20) ω − (cid:15) HF µ + iη (cid:21) d bµj , (56)where η is a positive constant of the order of the levelspacing of the µ -states. Since E ( t ) oscillates at the fre-quency ω P the ionization self-energy is dominated bythose terms in σ ( t − t (cid:48) ) that oscillate at energy (cid:15) HF µ (cid:39) ω P .We do therefore implement a frequency-independent ap-proximation ˜ σ abij ( ω ) ≈ ˜ σ abij ( ω P ), which in real time im-plies σ abij ( t − t (cid:48) ) = ˜ σ abij ( ω P ) δ ( t − t (cid:48) ). Substituting thisresult into Eq. (54) we getΣ > ion ( t, t (cid:48) ) = − iδ ( t − t (cid:48) )Γ ion ( t ) , (57)where Γ ion ,ij ( t ) = i (cid:88) ab E a ( t )˜ σ abij ( ω P ) E b ( t ) (58)is a self-adjoint positive-definite matrix for all times t .Thus, the approximate Σ ion is a local function of time asindicated in the left blue box of Fig. 5.The transient photoabsorption spectrum of the Kr gasin Fig. 4 [38] has been calculated using the 66 STO fromRef. [65] as basis, generating the input in System and
Ionization with the SMILES package [66, 67], freezing allelectrons below the 3 d shell and constructing Σ ion withthe HF states of positive energy. B. Kohn-Sham basis
In general the finite system of interest can be describedin terms of a single-particle basis formed by core statesand a rest. The rest is a set of active states, i.e., stateswith a population different from 0 or 1 because of dy-namical correlations or thermal fluctuations or externalfields. Let C and A be the set of core states and activestates respectively. Since, by definition, for i ∈ C ev-ery physically relevant many-body state is an eigenstateof ˆ c † i ˆ c i with eigenvalue 1, we can work in the truncated Hilbert space of many-body states having the core statesentirely filled. In this truncated Hilbert space the densitymatrix ρ satisfies again Eq. (43) but with a different HFHamiltonian h HF .To determine the HF Hamiltonian for the “active” elec-trons let us split the contributions (core and active) tothe Hartree and exchange potential V S H ,ij [ ρ ] ≡ (cid:88) mn ∈ S v imnj ρ nm , ij ∈ A (59) V S x ,ij [ ρ ] ≡ − (cid:88) mn ∈ S v imjn ρ nm , ij ∈ A . (60)where S = C , A and the indices i, j run in the active set A . Taking into account that ρ nm = δ nm for n, m ∈ C theequilibrium HF Hamiltonian in Eq. (15) can be rewrittenas h HF [ ρ ] = h eq+fc + V A H [ ρ ] + V A x [ ρ ] , (61)where h eq+fc = h eq + V C H + V C x (62)is the one-particle Hamiltonian plus the HF potentialgenerated by the frozen core electrons.So far we have not yet specified the single-particle ba-sis. We here consider the case of a Kohn-Sham (KS)basis. Hence we assume that electrons in the KS coreorbitals remain frozen and do not participate to the dy-namics. The equilibrium KS one-particle density matrixin the KS basis reads ρ KS ,nm = δ nm and the correspond-ing equilibrium KS Hamiltonian is diagonal and reads h KS = h eq + V C H + V xc + V A H [ ρ KS ] , (63)where V xc is the exchange-correlation potential of DensityFunctional Theory (DFT). In general, V C H + V xc is givenby the sum of the pseudopotential and the xc potentialgenerated by the active electrons. A comparison withEq. (61) allows us to express h eq+fc in terms of the KSHamiltonian according to h eq+fc = h KS − V xc − V A H [ ρ KS ] + V C x . (64)Depending on the system and laser pulse propertiesthe electrons in states with (cid:15) KS i < Λ max are explicitelypropagated through ρ , whereas states i = µ with en-ergy (cid:15) KS µ > Λ max are either assigned to S ion or S Auger ,see middle green box in Fig. 5 (here Λ max is a conver-gence parameter). Thus, CHEERS needs the KS eigen-values (cid:15) KS i (needed to construct the KS Hamiltonian h KS ,ij = δ ij (cid:15) KS i ), the matrix elements V xc ,ij , d ij and theCoulomb integrals v ijmn (needed to evaluate the Hartreepotential V A H [ ρ KS ] generated by the active KS electronsas well as the functionals V A H , V A x and Σ), the dipolematrix elements d iµ and KS energies (cid:15) KS µ with µ ∈ S ion (needed to calculate the ionization self-energy Σ ion ) and1the Coulomb integrals v Aijmµ and KS energies (cid:15) KS µ with µ ∈ S Auger (needed to calculate the Auger self-energyΣ
Auger and the kernel K ), see middle blue box in Fig. 5.This input contains the necessary quantities to constructthe functionals I tot and J , see Eqs. (43), as well as theHF Hamiltonian in Eq. (64). In fact, the only remainingunknown is V C x which, however, is usually small and canbe neglected. Of course, a non-negligible V C x does notintroduce extra complications for the CHEERS simula-tions. One could estimate this quantity by performingan all-electron KS calculation without pseudopotentials.The snapshots of the density variation of the pheny-lalanine aminoacid in Fig. 3 [40] has been calculatedby performing a DFT calculation with the QuantumEspresso package [68] using norm-conserving Troullier-Martins pseudopotentials [69] and the PBE approxima-tion [70] for V xc . The resulting KS states have then beenused to extract the matrix elements V xc ,ij , d ij , d iµ andthe Coulomb integrals v ijmn using the Yambo code [71](no Auger scattering was included since only valence elec-trons are ionized by the XUV pulse). C. An application to Argon: STO versus KS basis
As long as the single-particle basis is complete theCHEERS results are independent of the basis. The pur-pose of this Section is to illustrate this fact with an ex-ample. We consider the Argon atom and monitor thetime-evolution of the occupations of the HF orbitals af-ter a sudden ionization. The calculations are performedin two different basis:(i) the STO basis of Clementi-Roetti [72] consisting of 32basis functions. The input has been obtained with theSMILES package [66, 67] and no electron has been frozen(hence electrons of the K and L shells participate to thedynamics).(ii) the KS basis obtained by performing a DFT calcu-lation with the Octopus code [73] using norm-conservingTroullier-Martins pseudopotentials [69] and the Perdew-Zunger xc functional [74].In both cases we start from an initial density matrixcorresponding to the state of the system just after anionizing laser pulse has passed through the atom. Typi-cal attosecond pulses remove less than 1% of charge fromthe neutral system. Here, in order to highlight the ef-fects of correlations, we consider an initial density ma-trix ρ aσ,bσ (cid:48) (0) = δ σσ (cid:48) ρ a,b (0) where ρ (0) = ρ eq − δρ and δρ in HF basis reads δρ s, s = 0 . δρ p x , p x = 0 . δρ s, p x = δρ p x , s = − .
1. This corresponds to remove0.2 electrons of spin up and down.In Fig. 6 we compare the results in the two differentbasis as obtained by running CHEERS in the HF and 2Bapproximation. In HF the time-evolution is dominatedby 3 s ↔ p x transitions and resembles the evolution ofa noninteracting two-level system, in agreement with thefact that the HF theory is a single-particle theory and thesystem is weakly correlated. On the contrary, the corre- p x p x s s time (a.u.) o r b i t a l o cc upa t i on all electron - STO basis LDA + pseudopotential - Grid basis
FIG. 6: Time-dependent evolution of the occupations of theequilibrium HF orbitals 3 s and 3 p x after a sudden ionization,as described in the main text. The calculations have beenperformed using the HF (second and fourth panels) and 2B(first and third panels) approximation. lated 2B evolution highlights the occurrence of scatter-ings involving 3 p y and 3 p z electrons. In fact, the initialdensity matrix describes a mixture of charge neutral Arand multiply ionized Ar n + with n = 1 , . . . ,
6. In theconsidered Hilbert space Ar + can only give rise to theoscillation corresponding to the transition 3 s ↔ p whileAr can only give rise to oscillations corresponding tothe transitions 3 s ↔ s p and 3 s p ↔ p . These aredegenerate in HF although in reality they should not. 2Bcorrectly removes the degeneracy giving rise to the ob-served beating. Doubly and multiply ionized Ar atomscontributes less since we have removed only 40% of anelectron. Aside from the physical interpretation of theresults, the figure clearly show that the outcomes stem-ming from using two different basis (and procedures) arein a fairly good agreement.2 D. Retarded propagator
In this section we discuss the possible choices of theretarded Green’s function. The exact equation of motionfor G R reads (cid:20) i ddt − h HF ( t ) (cid:21) G R ( t, t (cid:48) ) = δ ( t, t (cid:48) )+ (cid:90) d ¯ t Σ Rtot ( t, ¯ t ) G R (¯ t, t (cid:48) )(65)to be solved with boundary condition G R ( t, t + ) = − i .The lowest order (in the Coulomb integrals) approxima-tion for G R is obtained by setting Σ Rtot = 0. In this casewe get the HF propagator G R ( t, t (cid:48) ) = − iθ ( t − t (cid:48) ) T e − i (cid:82) tt (cid:48) d ¯ t h HF (¯ t ) . (66)In CHEERS all approximations to G R have the form ofEq. (66) where h HF is replaced by some quasi-particleHamiltonian h qp = h HF + ∆. For ∆ = 0 we recover theHF propagator. The advantage of approximations like inEq. (66) is that for small δtG R ( t + δt, t (cid:48) ) (cid:39) e − i h qp( t + δt )+ h qp( t )2 δt G R ( t, t (cid:48) ) (67)and hence the calculation of G R ( t, t (cid:48) ) for all t (cid:48) < t scaleslinearly in t . Consequently, the overall scaling remainsquadratic with the maximum propagation time.The presence of a continuum due to leads and/or pho-toelectron states can be partially taken into account byapproximatingΣ Remb / ion ( t, t (cid:48) ) (cid:39) − ( i/ δ ( t − t (cid:48) )Γ emb / ion ( t ) (68)where Γ ion is defined in Eq. (58) whereasΓ emb ,ij = − (cid:88) kα T i,αk Im (cid:20) − (cid:15) αk + iη (cid:21) T αk,j , (69)see Ref. [21]. Setting Σ tot = Σ Remb + Σ
Rion in Eq. (65) onefinds ∆ = − ( i/ emb + Γ ion ).Correlation effects in the propagator can be taken intoaccount by making the approximation [21] (cid:90) d ¯ t Σ R ( t, ¯ t ) G R (¯ t, t (cid:48) ) (cid:39) (cid:20)(cid:90) d ¯ t Σ R ( t, ¯ t ) (cid:21) G R ( t, t (cid:48) ) ≡ ˜Σ( t ) G R ( t, t (cid:48) ) , (70)which amounts to add ˜Σ to h HF . We evaluate ˜Σ( t ) inEq. (70) using the GKBA and the adiabatic propagator˜ G Rad ( t, t (cid:48) ) = (cid:90) dω π e − iω ( t − t (cid:48) ) ω − h qp ( t ) + iη . (71)In this way we generate a self-consistent equation for˜Σ( t ) = ˜Σ[ ρ ( t ) , h qp ( t )]. In practice at the n -th timestep we determine ρ ( t n +1 ) from Eq. (43), then calcu-late h qp ( t n +1 ) using ˜Σ( t n +1 ) = ˜Σ( t n ), hence ˜ G R ( t n +1 , t (cid:48) )and finally the new ˜Σ( t n +1 ). This procedure is repeated a few times to achieve convergence. We point out thatpropagator used in the evaluation of I tot is G R and notthe adiabatic ˜ G R . The latter is only an auxiliary quantityto evaluate ˜Σ( t ).In CHEERS the quasi-particle Hamiltonian used for G R reads h qp = h HF − ( i/ α emb Γ emb + α ion Γ ion ) + α ad ˜Σ , (72)where the parameters α emb , α ion , α ad can be set to either0 or 1. VI. CONCLUDING REMARKS
To the best of our knowledge, CHEERS is cur-rently the only code which combines ab initio meth-ods with NEGF to calculate the nonequilibrium dynam-ics of molecular systems. CHEERS has already beenused to study the charge dynamics of molecular junc-tions [21], time-resolved photoabsorption of noble gasatoms [38, 75], charge separation in donor-acceptor com-plexes [39], charge migration in organic molecules [40]and time-resolved Auger decays [41]. The code handlesinputs obtained in any basis and can perform all-electronsas well as pseudopotential calculations.Currently, dynamical correlations are included at thelevel of the 2B approximation for the self-energy, al-though the scaling of the computational cost remainsidentical using a statically screened electron-electroninteraction for the exchange and second-order dia-grams [39]. Self-energy approximations like GW or T-matrix would restore the cubic KBE scaling with thenumber of time steps unless a GKBA for W or T is pro-vided, an advance which would be of utmost theoreticaland computational value.So far CHEERS simulations involving photoionizationprocesses have been performed using the KS contin-uum states generated either by the Quantum Espressocode [68] (planewave basis) or by the Octopus code [73](grid basis). Another promising alternative consists inusing a B-spline basis [76, 77]. In a recent work, theB-spline basis has been combined with the algebraic dia-grammatic construction method to calculate the attosec-ond pump-probe spectrum of carbon dioxide [78].Another limitation of CHEERS is that the nuclear po-sitions are kept fixed during the time evolution. Work toinclude harmonic effects through a Fan self-energy [79]evaluated with ab-initio frequencies and electron-nuclearcouplings is in progress. We are also planning to imple-ment the semiclassical Ehrenfest dynamics which requiresto calculate all one- and two-electron integrals along thenuclear trajectory. This extension of CHEERS is es-pecially relevant to access the subpicosecond timescale(10 ÷
100 fs) typical of charge transfer and charge separa-tion processes.3
Acknowledgements
We acknowledge inspiring and insightul discussionswith Emil B¨ostrom, Fabio Covito, Daniel Karlsson, Si-mone Latini, Andrea Marini, Yaroslav Pavlyukh, AngelRubio, Davide Sangalli, Anna-Maija Uimonen, Robertvan Leeuwen and Claudio Verdozzi. We also acknowl-edge EC funding through the RISE Co-ExAN (Grant No. GA644076). E.P. also acknowledges funding fromthe European Union project MaX Materials design at theeXascale H2020-EINFRA-2015-1, Grant Agreement No.676598 and Nanoscience Foundries and Fine Analysis-Europe H2020-INFRAIA-2014-2015, Grant AgreementNo. 654360. G.S. also acknowledge Tor Vergata Univer-sity for financial support through the Mission Sustain-ability Project 2DUTOPI. [1] F. Krausz and M. Ivanov, Rev. Mod. Phys. ,163 (2009), URL https://link.aps.org/doi/10.1103/RevModPhys.81.163 .[2] L. Gallmann, C. Cirelli, and U. Keller, Annual Review ofPhysical Chemistry , 447 (2012), URL https://doi.org/10.1146/annurev-physchem-032511-143702 .[3] M. Nisoli, P. Decleva, F. Calegari, A. Palacios, andF. Mart´ın, Chemical Reviews , 10760 (2017),pMID: 28488433, URL http://dx.doi.org/10.1021/acs.chemrev.6b00453 .[4] P. Danielewicz, Annals of Physics , 239 (1984).[5] G. Stefanucci and R. van Leeuwen, NonequilibriumMany-Body Theory of Quantum Systems: A ModernIntroduction (Cambridge University Press, Cambridge,2013).[6] K. Balzer and M. Bonitz,
Nonequilibrium Green’s Func-tions Approach to Inhomogeneous Systems (Springer,2012).[7] O. Konstantinov and V. Perel, SOVIET PHYSICSJETP-USSR , 142 (1961).[8] L. P. Kadanoff and G. A. Baym, Quantum statisticalmechanics: Green’s function methods in equilibrium andnonequilibirum problems (Benjamin, 1962).[9] L. V. Keldysh et al., Sov. Phys. JETP , 1018 (1965).[10] N.-H. Kwong and M. Bonitz, Phys. Rev. Lett. , 1768 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.84.1768 .[11] N. E. Dahlen and R. van Leeuwen, Phys. Rev. Lett. , 153004 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.98.153004 .[12] K. Balzer, M. Bonitz, R. van Leeuwen, A. Stan, and N. E.Dahlen, Phys. Rev. B , 245306 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.79.245306 .[13] P. My¨oh¨anen, A. Stan, G. Stefanucci, and R. vanLeeuwen, EPL (Europhysics Letters) , 67001(2008), URL http://stacks.iop.org/0295-5075/84/i=6/a=67001 .[14] P. My¨oh¨anen, A. Stan, G. Stefanucci, and R. vanLeeuwen, Phys. Rev. B , 115107 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.80.115107 .[15] M. P. von Friesen, C. Verdozzi, and C.-O. Almbladh,Phys. Rev. Lett. , 176404 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.176404 .[16] M. Puig von Friesen, C. Verdozzi, and C.-O. Almbladh,Phys. Rev. B , 155108 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.82.155108 .[17] K. Balzer, S. Bauch, and M. Bonitz, Phys. Rev. A , 022510 (2010), URL https://link.aps.org/doi/10.1103/PhysRevA.81.022510 .[18] K. Balzer, S. Bauch, and M. Bonitz, Phys. Rev. A , 033427 (2010), URL https://link.aps.org/doi/ 10.1103/PhysRevA.82.033427 .[19] M. Sch¨uler, J. Berakdar, and Y. Pavlyukh, Phys. Rev.B , 054303 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.054303 .[20] P. Lipavsk´y, V. ˇSpiˇcka, and B. Velick´y, Phys. Rev. B , 6933 (1986), URL https://link.aps.org/doi/10.1103/PhysRevB.34.6933 .[21] S. Latini, E. Perfetto, A.-M. Uimonen, R. vanLeeuwen, and G. Stefanucci, Phys. Rev. B ,075306 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.075306 .[22] K. Pernal, O. Gritsenko, and E. J. Baerends, Phys. Rev.A , 012506 (2007), URL https://link.aps.org/doi/10.1103/PhysRevA.75.012506 .[23] K. J. H. Giesbertz, O. V. Gritsenko, and E. J. Baerends,The Journal of Chemical Physics , 174119 (2010),URL https://doi.org/10.1063/1.3499601 .[24] K. J. H. Giesbertz, O. V. Gritsenko, and E. J. Baerends,The Journal of Chemical Physics , 18A517 (2014),URL https://doi.org/10.1063/1.4867000 .[25] M. Brics, J. Rapp, and D. Bauer, Phys. Rev. A ,013404 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.93.013404 .[26] M. Brics, J. Rapp, and D. Bauer, Journal of PhysicsB: Atomic, Molecular and Optical Physics , 144003(2017), URL http://stacks.iop.org/0953-4075/50/i=14/a=144003 .[27] F. Lackner, I. Bˇrezinov´a, T. Sato, K. L. Ishikawa,and J. Burgd¨orfer, Phys. Rev. A , 023412 (2015),URL https://link.aps.org/doi/10.1103/PhysRevA.91.023412 .[28] F. Lackner, I. B?ezinov, T. Sato, K. L. Ishikawa,and J. Burgdrfer, Journal of Physics: Conference Se-ries , 112084 (2015), URL http://stacks.iop.org/1742-6596/635/i=11/a=112084 .[29] S. Hermanns, N. Schl¨unzen, and M. Bonitz, Phys. Rev.B , 125111 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.125111 .[30] N. Schl¨unzen and M. Bonitz, Contrib. Plasma Phys. , 5 (2016), URL http://dx.doi.org/10.1002/ctpp.201610003 .[31] Y. Bar Lev and D. R. Reichman, Phys. Rev. B ,220201 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.220201 .[32] G. Pal, Y. Pavlyukh, W. H¨ubner, and H. C. Schneider,The European Physical Journal B , 327 (2011), URL https://doi.org/10.1140/epjb/e2010-10033-1 .[33] E. Perfetto, D. Sangalli, A. Marini, and G. Stefanucci,Phys. Rev. B , 205304 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.205304 .[34] D. Sangalli, S. Dal Conte, C. Manzoni, G. Cerullo, and A. Marini, Phys. Rev. B , 195205 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.195205 .[35] E. A. A. Pogna, M. Marsili, D. De Fazio, S. Dal Conte,C. Manzoni, D. Sangalli, D. Yoon, A. Lombardo, A. C.Ferrari, A. Marini, et al., ACS Nano , 1182 (2016),URL https://doi.org/10.1021/acsnano.5b06488 .[36] Sangalli, D. and Marini, A., EPL , 47004(2015), URL https://doi.org/10.1209/0295-5075/110/47004 .[37] E. Perfetto, D. Sangalli, A. Marini, and G. Stefanucci,Phys. Rev. B , 245303 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.245303 .[38] E. Perfetto, A.-M. Uimonen, R. van Leeuwen, and G. Ste-fanucci, Phys. Rev. A , 033419 (2015), URL https://link.aps.org/doi/10.1103/PhysRevA.92.033419 .[39] E. V. Bostr¨om, A. Mikkelsen, C. Verdozzi, E. Perfetto,and G. Stefanucci, Nano Lett. , 785 (2018), URL https://doi.org/10.1021/acs.nanolett.7b03995 .[40] E. Perfetto, D. Sangalli, A. Marini, and G. Stefanucci,The Journal of Physical Chemistry Letters , 1353(2018), URL https://doi.org/10.1021/acs.jpclett.8b00025 .[41] F. Covito, E. Perfetto, A. Rubio, and G. Stefanucci,Phys. Rev. A , 061401 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.061401 .[42] J. M. Luttinger, Phys. Rev. , A1505 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.135.A1505 .[43] F. G. Eich, A. Principi, M. Di Ventra, and G. Vignale,Phys. Rev. B , 115116 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.115116 .[44] F. G. Eich, M. Di Ventra, and G. Vignale, Phys. Rev.B , 134309 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.134309 .[45] F. Covito, F. G. Eich, R. Tuovinen, M. A. Sentef, andA. Rubio, Journal of Chemical Theory and Computation , 2495 (2018), URL https://doi.org/10.1021/acs.jctc.8b00077 .[46] M. Galperin, Chem. Soc. Rev. , 4000 (2017), URL http://dx.doi.org/10.1039/C7CS00067G .[47] Y. Meir and N. S. Wingreen, Phys. Rev. Lett. , 2512 (1992), URL https://link.aps.org/doi/10.1103/PhysRevLett.68.2512 .[48] A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B , 5528 (1994), URL https://link.aps.org/doi/10.1103/PhysRevB.50.5528 .[49] M. Sch¨uler and Y. Pavlyukh, Phys. Rev. B ,115164 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.115164 .[50] N. E. Dahlen and R. van Leeuwen, The Journal of Chem-ical Physics , 164102 (2005), URL https://doi.org/10.1063/1.1884965 .[51] K. Balzer, S. Hermanns, and M. Bonitz, EPL (Euro-physics Letters) , 67002 (2012), URL http://stacks.iop.org/0295-5075/98/i=6/a=67002 .[52] N. S¨akkinen, M. Manninen, and R. van Leeuwen, NewJournal of Physics , 013032 (2012), URL http://stacks.iop.org/1367-2630/14/i=1/a=013032 .[53] M. Hopjan, D. Karlsson, S. Ydman, C. Ver-dozzi, and C.-O. Almbladh, Phys. Rev. Lett. ,236402 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.236402 .[54] Y. B. Lev and D. R. Reichman, EPL (Europhysics Let-ters) , 46001 (2016), URL http://stacks.iop.org/ 0295-5075/113/i=4/a=46001 .[55] N. Schl¨unzen, J.-P. Joost, F. Heidrich-Meisner, andM. Bonitz, Phys. Rev. B , 165139 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.165139 .[56] A.-M. Uimonen, E. Khosravi, A. Stan, G. Stefanucci,S. Kurth, R. van Leeuwen, and E.K.U. Gross, Phys. Rev.B , 115103 (2011).[57] C.-O. Almbladh, A. L. Morales, and G. Grossmann,Phys. Rev. B , 3489 (1989), URL https://link.aps.org/doi/10.1103/PhysRevB.39.3489 .[58] C.-O. Almbladh and A. L. Morales, Phys. Rev. B , 3503 (1989), URL https://link.aps.org/doi/10.1103/PhysRevB.39.3503 .[59] M. Cini, Solid State Communications , 681(1977), URL .[60] G. A. Sawatzky, Phys. Rev. Lett. , 504 (1977), URL https://link.aps.org/doi/10.1103/PhysRevLett.39.504 .[61] G. Stefanucci, S. Kurth, A. Rubio, and E. K. U. Gross,Phys. Rev. B , 075339 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.075339 .[62] G. Stefanucci, E. Perfetto, and M. Cini, Phys. Rev. B , 115446 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.81.115446 .[63] A. Kalvov´a, B. Velick´y, and V. ˇSpiˇcka, EPL (EurophysicsLetters) , 67002 (2018), URL http://stacks.iop.org/0295-5075/121/i=6/a=67002 .[64] E. Goulielmakis, Z.-H. Loh, A. Wirth, R. Santra,N. Rohringer, V. S. Yakovlev, S. Zherebtsov, T. Pfeifer,A. M. Azzeer, M. F. Kling, et al., Nature , 739 (2010).[65] C. Bunge, J. Barrientos, and A. Bunge, Atomic Dataand Nuclear Data Tables , 113 (1993), ISSN 0092-640X, URL .[66] J. Fern´andez Rico, I. Ema, R. L´opez, G. Ram´ırez and K.Ishida, in Recent Advances in Computational Chemistry:Molecular Integrals over Slater Orbitals , eds. T. Ozdoganand M. B. Ruiz (Transworld Research Network, 2008),pp. 145.[67] J. F. Rico, R. L´opez, I. Ema, and G. Ramrez, Jour-nal of Computational Chemistry , 1987 (2004),URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.20131 .[68] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, et al., Journal of Physics: Condensed Mat-ter , 395502 (2009), URL http://stacks.iop.org/0953-8984/21/i=39/a=395502 .[69] N. Troullier and J. L. Martins, Phys. Rev. B , 1993 (1991), URL https://link.aps.org/doi/10.1103/PhysRevB.43.1993 .[70] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865 .[71] A. Marini, C. Hogan, M. Gr¨uning, and D. Varsano, Com-puter Physics Communications , 1392 (2009), URL .[72] E. Clementi and C. Roetti, Atomic Data andNuclear Data Tables , 177 (1974), URL .[73] X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas,I. Theophilou, N. Helbig, M. J. Verstraete, et al., Phys.Chem. Chem. Phys. , 31371 (2015), URL http://dx.doi.org/10.1039/C5CP00351B .[74] J. P. Perdew and A. Zunger, Phys. Rev. B ,5048 (1981), URL https://link.aps.org/doi/10.1103/PhysRevB.23.5048 .[75] E. Perfetto, A.-M. Uimonen, R. van Leeuwen, andG. Stefanucci, Journal of Physics: Conference Se-ries , 012004 (2016), URL http://stacks.iop.org/1742-6596/696/i=1/a=012004 .[76] M. Ruberti, V. Averbukh, and P. Decleva, The Journal of Chemical Physics , 164126 (2014), URL https://doi.org/10.1063/1.4900444 .[77] M. Ruberti, P. Decleva, and V. Averbukh, Phys. Chem.Chem. Phys. , 8311 (2018), URL http://dx.doi.org/10.1039/C7CP07849H .[78] M. Ruberti, P. Decleva, and V. Averbukh, Journal ofChemical Theory and Computation , null (0), URL https://doi.org/10.1021/acs.jctc.8b00479 .[79] H. Y. Fan, Phys. Rev. , 900 (1951), URL https://link.aps.org/doi/10.1103/PhysRev.82.900https://link.aps.org/doi/10.1103/PhysRev.82.900