A correlated-polaron electronic propagator: open electronic dynamics beyond the Born-Oppenheimer approximation

In this work we develop a theory of correlated many-electron dynamics dressed by the presence of a finite-temperature harmonic bath. The theory is based on the ab-initio Hamiltonian, and thus well-defined apart from any phenomenological choice of collective basis states or electronic coupling model. The equation-of-motion includes some bath effects non-perturbatively, and can be used to simulate line- shapes beyond the Markovian approximation and open electronic dynamics which are subjects of renewed recent interest. Energy conversion and transport depend critically on the ratio of electron-electron coupling to bath-electron coupling, which is a fitted parameter if a phenomenological basis of many-electron states is used to develop an electronic equation of motion. Since the present work doesn't appeal to any such basis, it avoids this ambiguity. The new theory produces a level of detail beyond the adiabatic Born-Oppenheimer states, but with cost scaling like the Born-Oppenheimer approach. While developing this model we have also applied the time-convolutionless perturbation theory to correlated molecular excitations for the first time. Resonant response properties are given by the formalism without phenomenological parameters. Example propagations with a developmental code are given demonstrating the treatment of electron-correlation in absorption spectra, vibronic structure, and decay in an open system.


I. INTRODUCTION
The small-polaron transformation of the electronic Hamiltonian was originally developed in the 1960's 91 , and more recently revived 92,93 in the many-electron context. It is a classic example 94 of the utility of canonical transformations in quantum physics. Its usefulness is well-established 94 yet it is also experiencing renewed interest [95][96][97][98] . In particular, second-order master equations in the polaron frame afford good results in all bath strength regimes [99][100][101] employing the variational technique of Harris and Silbey 102 . In the many-electron case, there has also been some recent pioneering work toward developing random-phase approximation equations 92,93 . The electronic structure community has also produced some related work [103][104][105][106][107][108] , including phenomenologically damped 109 response theory.
In this paper we formulate and implement a correlated many-electron master equation that overcomes several limitations of the adiabatic Born-Oppenheimer approximation, and includes effects such as excited-state lifetimes and vibronic structure. The basic goal of our method is to produce electronic spectra for small molecules that include the effects of coupling to an environment. With continued development our formalism could describe environmental localization of electronic states and the decay of quantum entanglement in correlated electronic systems.
The theory we present exploits the polaron transformation to combine both electron correlation and system-bath couplings in a single perturbation theory. In the transformed frame, high-rank quantum expressions are dressed by environmental factors, which cause them to decay during dynamics. This introduces the possibility for environmentally induced decay of the correlations in an electronic system, thereby making the problem computationally more tractable. We also discuss how a general one-particle perturbation to the electronic system may be treated in a closely related untransformed version of the theory. Physically relevant coupling models of this form are numerous, and several examples include nuclear motion coupled to electronic degrees of freedom 110,111 , Coulomb coupling to a nearby nanoparticle surface 112 , the electromagnetic vacuum 113 , and perhaps even Coulomb coupling to surrounding molecules in a condensed phase.
The present method is distinguished from previous work by a few characteristic features.
Unlike virtually all master equation approaches, it treats the dynamics without assuming the many-body problem of the electronic system to be solved. We emphasize that our correlated theory is one half of a final theory with the other involving a recipe to calculate system specific couplings with Ehrenfest type dynamics or another scheme. The reason being that, in general, the appropriate basis of a few electronic states to prepare a master equation can't be found a priori. As a result, this theory begins from a system-bath Hamiltonian which is well-defined and atomistic in terms of single-electron states and employs a non-Markovian equation of motion (EOM) in place of phenomenological damping.
Conveniently, we find that high-rank operator expressions responsible for the computational intractability of exact, closed, many-particle quantum mechanics are multiplied by factors which exponentially vanish in many circumstances.
The picture of electronic dynamics offered by a master equation is complementary to the time-dependent wave-packet picture of absorption [114][115][116][117][118][119][120] . In the latter approach the electronic degrees of freedom are described as a superposition of a few adiabatic surfaces, and nuclear wave-packets are studied in the dynamics. In our approach the electronic degrees of freedom are described in all-electron detail, but the dynamics of the nuclei are approximated by the choice of the Holstein Hamiltonian. It is possible to combine a wave-packet calculation of the bath correlation function with the untransformed equation of motion presented in this work, to leverage the physical strengths of both approaches.
We demonstrate the implementation of the formalism in a pilot code, and apply it to calculate some dynamic properties of model small molecules. The spectra produced by the electron correlation theory are shown to compare favorably to related methods in the adiabatic limit. Vibronically progressed spectra are shown to be produced by the dressed theory, and a Markovian version is applied to a basic model of electronic energy transport between chromophores. Finally the undressed, uncorrelated theory is used to simulate the ultraviolet absorption spectrum of 1,1-diflouroethylene and compared with available experimental data.

A. Hamiltonian
As in the work of references 92,93 , we use the non-relativistic ab-initio many-electron Hamiltonian with a Holstein-type 121 (linear) coupling to a bath of non-interacting bosons (summation over repeated indices is implied throughout this paper unless stated otherwise), given Here, a † s creates an electron in the single-particle basis state s, while b † k creates a bosonic bath particle in the kth mode.F is the Fock operator of the reference determinant, andV is the two-electron part of the electronic Hamiltonian in the same single-particle basis. The third term is the boson Hamiltonian,Ĥ b and the last is the coupling of the electronic system to the bath. For a general bath mode with dimensionless displacement Q k , the bi-linear Hamiltonian. In the following discussion, we will assume the one-electron parts ofĤ andF to be diagonal in the (canonical) one-electron basis with eigenvalues p = f p p . We now introduce the displacement operator, which generates the polaron transformation. Since the operatorŜ is anti-Hermitian, it generates a unitary transformation of the electron-boson HamiltonianĤ →H = e −ŜĤ eŜ.
Explicitly, the polaron transformed Hamiltonian is given bỹ where the one-particle part of the transformed Hamiltonian is given bỹ and we have introduced the reorganization energy λ p = k M p k /ω k . The transformed electron-electron interaction is given bỹ where the transformed matrix elements arẽ and we have defined the bath operatorsX p aŝ For future reference it will also be useful to define dressed electronic creation and annihilation operators, denoted byã † s ≡ a † s X † s andã s ≡ a s X s . The key feature of the polaron transformation is thatH has no electron-phonon coupling term. As a result, the two-electron and electron-boson parts of the original Hamiltonian are combined into a single term which now couples two dressed electrons and two dressed bosons (Eq. 6). One should intuitively imagine the situation depicted in Fig. 1, where two displaced electronic energy surfaces are dragged into alignment by the polaron transformation, altering the coupling region between them, which now absorbs the electron-boson coupling.
In what follows, our expressions will be derived in the interaction picture with respect to Eq. 6 and then switched to the Schrödinger picture. The harmonic nature of the bosons means that the correlation functions ofX operators (in any combination and at multiple times) can be given as simple functions of ω p ,M and β = k b T 92 . The transformed electronboson problem takes the form of the usual many-electron problem with a time-dependent electron-electron interaction. This allows us to harness powerful methodologies that stem from two distinct areas of research: quantum master equations 123 and quantum chemistry methods 124 . We exploit this feature to produce a model of electronic dynamics which treats system-bath dynamics and correlation effects within the same perturbation theory.

B. Electronic dynamics
Having reviewed the polaron transformation in the previous section, we now combine it with the time-convolutionless perturbation theory (TCL) 125,126 to arrive at the central theoretical result of this manuscript. Our goal is to derive an equation-of-motion for a dressed particle-hole excitation operator, which we denoteõ i a ≡ o a iã † aã i (here i, j, k... are zeroth-order occupied levels and a, b, c... are unoccupied). We consider only the particle-hole part of the density matrix, commonly referred to as the Tamm-Dancoff approximation 127,128 (TDA), to simplify the derivation of our equations and avoid the possible issues associated with the non-linearity of the other blocks 129 . To derive Fock-space expressions 130 , it is convenient to assume the initial equilibrium state is the canonical Hartree-Fock (HF) determinant i.e. the initial density matrix is |Ψ(0) Ψ(0)| ≈ |0 0|, with |0 the HF determinant. This assumption relies on two approximations: 1) The first excited-state energy of the systems we will study is much larger than k b T , so that the system is effectively in the ground-state.
2) The initial state of the system is weakly correlated and hence dominated by a single Slater determinant. Assumption 1 is clearly valid for the small molecular systems we will study, while assumption 2 requires more care and the effects of initial correlations are discussed in detail in section II. E.
(9) P F is a Fock space projector onto maps between single-electron {a † a} operators like the typical projector of a partitioned electron-correlation perturbation theory ie: This partitioning is consistent with the perturbative ordering ofH by powers ofṼ . We treat the description of the many-electron state in terms of only one-body operators (neglecting all higher density matrices) on the same footing as tracing over the bath degrees of freedom 131 , so both phenomena are easily incorporated in the same master equation.
The effective Liouvillian becomes time-dependent due to the polaron transformation.
Consequently, there exists no simple analytical formula for its Fourier transform. Instead, we must give a differential representation of the EOM for a particle-hole excitation, which can then be integrated numerically. Given these projectors, the time-convolutionless perturbation theory 125 (TCL) over the space of P is: 132 .
This is written in the interaction picture where L = −i[Ṽ (t), ·]. The first term above is the uncorrelated part of the evolution 133 , the second is a homogenous term reflecting correlation between the system and the bath, and the last term is an inhomogeneity reflecting correlations of the initial state. The interaction picture perturbation 134 is . Expanding the first term over the one-particle space and moving into the Schrodinger picture we obtain (left arrows, ←, are used to indicate the contribution of a term to dõ/dt ): We have introduced a shorthand for the correlation function B Term (13)  The second homogeneous term introduces bath correlation functions between boson operators occurring at different times (t and s) according to Moving the electronic part into the Schodinger picture and rearranging this becomes: Expanding the commutators in Eq. 17, applying Wick's theorem to remove the many vanishing terms, and enforcing the connectivity constraint, one obtains many topologically distinct terms. In our implementation this is done automatically before the execution of a simulation. The terms are easily related to terms which occur in the expansion of the second-order Fermion propagator (SOPPA) [136][137][138] and diagonalization-based excited-state theories like CIS(D) 139 . Since we employ the time-convolutionless perurbation theory, the oscillating e i∆t factors are different than those which occur in Rayleigh-Schrodinger perturbation theories (in the energy-domain denominator 1 ω−∆ ) and warrants further study. We some terms from a hole → particle excitation which is obtained by multiplying both sides of equation (27) on FIG. 2: The particle-hole → particle-hole skeleton Brandow 141 diagrams in the second-order homogeneous term, the first two come from the V t V s o term in the expanded double-commutator and third from V t oV s . The fourth occurs in V t oV s , V s oV t and V t V s o. The boson correlation function depends on any lines emerging from the ellipses which represent virtual electron scattering events at times t, s.
This term has six indices overall, but can be factorized as follows: The calculation of I a c scales fifth order with the number of single electron states, and linearly with the number of bath modes (which go into the calculation of B(t)). This low scaling with number of bath modes is inherited from the approach of Silbey 102 , Jang 97 and Nazir 140 .
The algebraic version of the second term in Fig. 2 is: This term is sixth order with the size of the system, with the boson correlation function preventing a desirable factorization of the V-V contractions. 14 sixth-or-less-order terms are found which couple hole-particle excitations to each other; skeletons 141 of these are shown diagrammatically in Fig. (2) with explicit expressions for each of these terms given in the supplementary information 142 .We should note that our formalism lacks several terms which occur in the SOPPA because of the Q projector and the absence of V t o vac oV s ordered terms that are invoked in the formal expansion of the time-ordered exponential.
As described so-far the EOM is not time-reversible even whenM = 0 because neither the homogenous term nor the normal excitation operator take the form of an anti-Hermitian matrix. 143 . We remedy this by adding the terms which result from the normal excitation's operator's Hermitian conjugate, ie: the indices of the vacuum and o are swapped, and the signs of ∆ s and prefactor are flipped and the new term is added to the other perturbative terms. For example Term 20 becomes the following two terms: Making this modification, the linear response spectrum of the adiabatic model is not significantly altered, but the adiabatic norm conservation is enforced, and changes by less than 1 percent after 1700 atomic units in the case of H 4 with 4 th order Runge-Kutta and a timestep of 0.05 au.

C. Dipole Correlation function
Like any canonical transformation 144 , the polaron transformation preserves the spectrum of the overall electron-phonon Hamiltonian but the statistical meaning of the state related to a particular eigenvalue is changed. In other words: the operators of our theory are different objects from the adiabatic Fock space. To lowest order in system-bath coupling 140 , electronic observables like the time-dependent dipole correlation function, which predicts the results linear optical experiments, can be obtained by generating the transformed property operator Within the Condon approximation the dipole correlation function is then the product of the electronic trace and bath trace over theX operators introduced inμ. This adds a bath correlation function to the usual observable. Taking all the relevant expectation values gives the following explicit numerical formula for the dipole moment expression: This expression for the dipole-dipole correlation function assumes the bath remains at equilibrium. CIS likewise only offers a zeroth-order oscillator strength 145 . The spectra in this work are generated by "kicking" the electronic system with the dipole operator instantaneously, and Fourier transforming the resulting dipole-dipole oscillations.

D. Untransformed Version
Because of the polaron transformation the above formalism is accurate in the strong bath regime with diagonal system-bath couplings. To treat an off-diagonal, weak coupling one can develop the complementary untransformed theory. With a bi-linear system-bath coupling of the form: development of an uncorrelated particle-hole equation of motion follow similarly to the above with H sb taking the place ofṼ . The untransformed version also makes it possible to introduce an Ehrenfest scheme for the nuclear bath which may be pursued in future work.
The projector of the untransformed version is simply the equilibrium trace over b operators.
The TCL produces second-order contribution of the form: is given by the usual formula 147 : Furthermore, we find that the reorganization energy is given by λ i = ω i d 2 ii which can be subtracted from the original Hamiltonian.

E. Initial Correlations
In the above we have expanded on existing master equations by giving the electronic system a many-body Hamiltonian with Fermionic statistics. We have also improved typical electronic response theories by incorporating bath dynamics, but we have assumed that the initial state was a single determinant and employed the corresponding Wick's theorem.
These approximations are made in many other treatments of electronic spectra 110,148,149 , and are acceptable for situations where a molecule begins in a nearly determinantal state.
Optical absorption experiments of gapped small molecules belong to this regime. To rigorously study electron transport in a biased junction 150 , or otherwise more exotic initial state requires a treatment of initial correlations [151][152][153][154] , as in the non-equilibrium Green's function(NEGF) method [155][156][157] . These methods first propagate an initial determinant in imaginary time 150 thermalizing and correlating the system before the dynamics, but numerical applications of the NEGF formalism require storage an manipulation of a state variable with three indices: G pq (ω), and are usually limited to small systems 158 , and short times. We can similarly perform an imaginary time propagation, and have performed some exploratory calculations which are described in the supporting supporting material 142 .
Instead of performing an imaginary time integration, one can replace the usual Wick's theorem, and build the theory beginning from a state which is already correlated. Extended normal ordering 159 , makes it possible to take expectation values with a multi-configurational reference function via a generalization of Wick's theorem 160 . Using the extended normal ordering technique, the expressions of the present paper can be promoted to treat a general initial state directly, so long as the density operator of that state is known. This approach would eliminate the need for any adiabatic preparation of the initial state, but would generate more complex equations, and is a matter for future work.
For weakly correlated systems it's reasonable to keep the approximation that the usual Wick's theorem applies to the initial state of the electronic system, and instead of adiabatically preparing an initial state, choose a perturbative approximation to the correlated part of the initial state, Qõ(0) = 0 −∞ dsL(s)õ(0), and treat the first inhomogenous 161 term of the master equation: The bath parts of the inhomogenous term are known to be relatively unimportant from studies of tight-binding Hamiltonians 97 , and only slightly perturb the results of a propagation for short times, and so we will not include them.
It is interesting to note that with the addition of this inhomogenous term, there is a (ω), with the dimension of the particle-hole space: In the aboveT 2 is the second-order excitation operator from Moller-Plesset perturbation theory. These terms correspond respectively to the Fourier transform of our PVP, PVQVP, and I(t) terms with different denominators. This correspondence suggests the present method should produce linear response spectra of quality similar to CIS(D), which is usually slightly Approximation Extension |Ψ eq ≈ |0 Extended normal ordering 159

Second-Order inṼ TCL-4 162 µ(t) ≈ T r(μ ·õ(t))
Use basis commuting with µ TDA Include the additional blocks 145 . Orbital relaxation Variational condition 102   TABLE I: Limitations of this work and how they may be relaxed. Orbital relaxation isn't really an approximation, per-se, but would be beneficial for the results of the perturbation theory.
better than TDDFT. Because the formalism is somewhat involved and many approxima-

C. Energy Transport and Markovian Evolution
The continuous integration of rank-6 bath correlation tensors required for the Non-Markovian propagation executed above is a rather costly proposition for large systems. This is especially true if one would like to study incoherent electronic energy transport which takes place in times on the order of picoseconds, roughly a million times more than the electronic timestep required to integrate Eq. (27) for a typical molecule. A useful approximation to overcome this is a Markov approximation, by which we mean taking the limit of the integral to infinity for each term, such as: R ab,cj ij,ak = lim t→∞ t t 0 B ab,cj ij,ak (t, s)(e i(∆ cj ab )(t−s) )ds (29) where R is now a factor replacing the integral expression. Each term possesses its own R tensor, but these must only be calculated once. Numerically, this limit can be taken in the case of a continuous super-ohmic bath with cutoff frequency ω c by introducing a cutoff time t c , above which the bath correlation function is assumed to be equal to it's equilibrium value, which is a good approximation for βω c (t c ) 2 >> 2.
We have used the letter R to suggest the analogy between this time-independent rate tensor and that occurring in the Redfield theory 169,170 , although in this theory there are 28 such fifth and sixth rank tensors. These tensors also differ because they use the transformed,  The evolution of this "coherence" isn't captured well in the Tamm-Dancoff approximation because we have truncated the blocks of the propagator coupling particle-hole excitations and their conjugate hole-particle modes.

D. Untransformed and uncorrelated version
The transformed version of the theory has several interesting formal advantages, but to treat larger systems or off-diagonal system-bath couplings we have also presented the Accurate protocols for providing a per-orbital spectral density must be developed and the accuracy of the resulting lineshapes and lifetimes must be assessed. Simply projecting the orbital energy gradient onto the normal modes of the molecule is an obvious choice, but it would be appealing to treat the influence of surrounding molecules in a similar way.
The Tamm-Dancoff equation of motion provided in this work propagates only one-block of a much larger transition density matrix which must be treated to capture coherence phenomena between particles amongst themselves (and holes). The treatment of a thermalized initial condition is an interesting, but computationally demanding question. The payoff for developing these additional features would be electronic excited states with more realism than can be offered in the adiabatic picture. These states would be naturally localized, with a size depending on vibronic structure and temperature and evolve and relax irreversibly.

V. APPENDIX: EXPRESSIONS FOR THE CORRELATION FUNCTIONS AND FACTORIZATION
Time-ordered harmonic correlation functions (HCF's) of all orders of the type XX † ...XX † were made available in Dahnovsky's pioneering work 92 by expanding the exponentials inX as power series inÂ k = (b k e iω k t −b † k e −iω k t ) and applying Wick's theorem to the resultingÂ operator strings paying careful attention to the combinatorial statistics. Since we employ a master equation theory rather than a Green's function theory, our version of the HCF depends on the sign of (t − s) but is otherwise the same after making the simplifications which appear in our second order theory.
If ∆ < 0 then relaxation occurs ifM at tMat s > 0. In most applications of Master equations a single, positive, spectral density is assumed for all states which basically parameterizes M at tMat s as a function of ω. This approximation is more than a mere convenience; ifM is assigned generally and different states are allowed to couple to a single frequency with different strengths it's quite easy to for some elements of the EOM to haveM at tMat s > 0 causing exponential growth in the Markovian limit and at short times. We note that the requirement that the correlation function be easily integrable is only a pre-requisite for the polaron transformation. If only the time-convolutionless equation of motion is used, any correlation function which is known can be easily incorporated.
Without further approximation a sixth-order number of B's must be calculated and integrated (in our code a third order Gaussian Quadrature with the electronic time-step was found sufficient), since at least two indices are shared between eachṼ . Consider the cost limiting term Eq. (20). It would be advantageous to evaluate the implied sum over c and make a 5-index intermediate, since then the remaining indices are also a fifth order loop, unfortunately the HCF depends on index c. On the other hand, High-rank operator strings are exponentially damped by the presence of this HCF in the Markovian limit. It seems likely that this feature could be used as a new locality principle which would lift the curse of dimensionality in the strong bath regime.