Time-Dependent Current-Density Functional Theory for Generalized Open Quantum Systems

bond


Introduction
A closed system is a quantum mechanical state that evolves under Hamiltonian evolution, therefore obeying Schro¨dinger's equation.In practice, however, a quantum system is not closed but interacts with the environment, and for many physical situations of interest, this interaction must be properly addressed.As an example, environmental effects are central in quantum decoherence and quantum thermodynamics, and the framework to study them is open quantum systems (OQS). 1,2With the increasing possibility of designing, manipulating and controlling objects at the nanoscale, manybody OQS are becoming ubiquitous subjects of study in current research spanning a broad range of disciplines in the physical sciences from traditional condensed matter 3 and chemical physics, 4 to the emerging fields of biophysics 5 and quantum information science. 6In order to achieve a substantial understanding of these systems, accurate yet computationally tractable theoretical techniques to study their time evolution are required.Over the last fifty years, significant progress in many-body theory has resulted from the introduction of density functional theory techniques both in the time-independent (DFT) 7,8 and dependent (TD-DFT) 9 domains.In particular, incorporating the OQS formalism into TD-DFT would provide a convenient set of tools for studying a vast number of dynamical processes such as excitations of molecules embedded in complex biological environments, 10 spin diffusion, 11 molecular conductance, 12 particle thermalization, 13 and many other interesting phenomena.In this article, we shall mostly concern ourselves with TD-DFT and its variants, and therefore only mention ground state DFT when needed.
Broadly speaking, TD-DFT reformulates time-dependent quantum mechanics in terms of particle densities instead of wavefunctions (or density matrices), thus allowing for more affordable computational scaling than standard many-body theories when it comes to the resources needed to study the time evolution of a closed system.Whereas the calculation of wavefunctions depends on 3N spatial variables (6N in the case of density matrices) and time, with N being the number of particles in the system, the calculation of particle densities depends only on three spatial and one time variables.However, we should emphasize that particle densities do not render the many-body problem trivial.In fact, the accurate reproduction of the original system's particle densities using auxiliary non-interacting systems, a strategy known as the Kohn-Sham (KS) scheme, still requires considerable work in the crafting of good functional approximations for the effective scalar KS scalar potentials via the development of the so-called exchange-correlation (xc) potentials.
So far, most of the development of TD-DFT has occurred in the context of closed systems.We are only aware of a few attempts to extend TD-DFT to the treatment of OQS.Gebauer, Car, and Burke (GCB) 14 have proved a Runge-Gross (RG) theorem 9 to include Markovian OQS of the Kossakowski-Lindblad form into TD-DFT.Di Ventra and D'Agosta (DADV) 16 have performed a similar adaptation for TD-CDFT in the context of stochastic Schrodinger equations.From a slightly different perspective, Chen and co-workers 15 have carried out a computational study of a dissipative molecular device using regular TD-DFT; they justified the validity of this formalism to treat OQS by invoking the holographic electron density theorem.Our present investigation is closely related to the first two efforts above.In particular, we will devote this article to TD-CDFT since our study on TD-DFT and OQS is carried out elsewhere. 17D-CDFT differs from TD-DFT in the fact that it is a theory based on current densities rather than particle densities.
Whereas TD-DFT is formally exact, its most common practice is based on spacially local functional approximations which have been successful in many applications (see ref. 18 and 19), but have also failed in some situations such as with the description of charge transfer excited states. 20On the other hand, TD-CDFT is a more expensive theory compared to TD-DFT due to the intrinsic vectorial nature of the current density (it requires 6 spacial variables plus time), 21,22 but it has been pointed out that the problems of ultranonlocality in space that pervade TD-DFT do not appear in this theory. 23ith this in mind, spacially local exchange correlation vector potentials 24 have been implemented within TD-CDFT to successfully predict polymer polarizabilities, 25 electronic properties of weakly disordered systems, 26 and recently, charge and spin dynamics in ultracold gases, 27 among many other applications.Our goal is to generalize the framework of TD-CDFT so that the set of problems consisting of either Markovian or non-Markovian OQS can be studied within the theory.We will do so by establishing formal results on the RG theorem and on the representability of open systems with auxiliary KS systems.
The remainder of this article is organized as follows: in the next section, we generalize the RG theorem for arbitrary OQS, comment on the implications it has with respect to TD-CDFT, and note on the differences between standard TD-CDFT and this OQS version.In section 3, we relate the construction given in section 2 with the problem of continuity and master equations, and thus conclude that the proof of the existence of the KS scheme suggested by DADV might have a limited applicability; however, we suggest a way to improve this problem in section 4. In section 5, we suggest a novel KS scheme where instead of using an open system as the auxiliary KS system, we work with a closed one.This constitutes the main result of our article, since we suggest not only switching from an interacting to a non-interacting system, but also from an open to a closed system.We comment on this strategy comparing it with the previously suggested schemes, and on issues related to the development of functionals for this particular proposal.Finally, in section 5, we give a summary of the article together with conclusions and comments about future work.

Runge-Gross theorem for generalized open quantum systems
The most standard way to treat OQS 1 is through the density matrix r(t) with the differential form of its evolution given by a master equation of the form: where the right hand side of the equation consists on a closed-system evolution by the Hamiltonian H ˆS(t), plus a dissipative evolution characterized by a memory kernel K(t,t 0 ).We note that although any open quantum system evolution can be written in the form (1), 28,29 not every memory kernel K(t,t 0 ) denotes a valid physical evolution, and in fact, additional constraints must be applied to it in order for it to have any meaning. 2,30Although the exact dynamics can still be recasted in this form, in practice, due to the impossibility in most practical cases to characterize the bath and its interaction with the system in an explicit way, one is resorted to simple models for K(t,t 0 ), thus rendering eqn (1) as a reasonable approximation to the real evolution of the system of interest.
The simplest class of master equations has the semi-group property under the Markov approximation; 31,32 this class assumes an ideal memoryless environment that does not act back on the system.The most common example of this situation is the Kossakowski-Lindblad (KL) form, [31][32][33] where the action of the memory kernel on the density matrix is given by describes the effects on the system of an ideal bath with L ˆi(t) + being a KL jump generator.Here a ij are the real valued jump rates and R is the dimension of the many-body Hilbert space of the system in consideration.
On the other hand, there are also several formalisms for treating non-Markovian equations in order to account for more realistic environments with memory.The projector operator method splits the evolution into the relevant part, the system, and the rest, the environment in a manner that makes it tractable to account for effects of the bath on the system. 28,29This method can be applied to derive the so-called time convolutionless master equations which incorporate memory effects that are local in time. 2 Their connection to non-Markovian dynamical maps was studied in ref. 30.In any case, in this investigation we only need to assume that K(t,t 0 ) will be of a form with a valid physical interpretation without subscribing to a particular class of master equations.Now, the unitary part of the evolution is characterized by the many-body Hamiltonian H ˆS(t) which contains both vector and scalar potentials, Ã(r,t) and Ṽ(r,t), respectively, and an interparticle pairwise symmetric potential U(r i ,r j ): Without loss of generality, we can set V(r,t) = 0 via the corresponding gauge transformation (see ref. 22), and we shall assume that this transformation has been performed hereafter.Also, just for nomenclature purposes, we define the particle density operator as nðrÞ ¼ P i dðr À ri Þ and the current density operator as ĵðr; tÞ ¼ indicates a trace with respect to r(t).We emphasize that the time dependence of the expectation values hO ˆ(r,t)i t will stem both from the explicit time dependence of the operator O ˆ(r,t) and from the evolution of r(t) due to the master eqn (1). 34ow, just as in the original Runge and Gross (RG) article 9 we shall define certain mappings in order to establish the language of what will be our TD-CDFT for OQS (see Fig. 1).Working with a fixed initial density matrix r(0), interparticle potential U(r i ,r j ), and memory kernel K(t,t 0 ): 1. We first establish a map from vector potentials to density matrices F: Ãðr; tÞ !rðtÞ by simply using the master equation (1).
2. Next, by computing the particle and current densities at a given point in space and time, we obtain a map between vector potentials and coordinate pairs of particle and current densities G: Ãðr; tÞ !ðhnðrÞi t ; h ĵðr; tÞi t Þ.In the next paragraph we will establish the OQS version of the RG theorem, namely, that G is one-to-one and thus invertible from the image set of G.This fact allows us to write rðtÞ ¼ FG À1 ðhnðr Þi t ; h ĵðr; tÞi t Þ, and therefore, establish that the expectation value of any observable hO ˆ(r,t)i t is a functional of r(0), U(r i ,r j ), K(t,t 0 ), hnˆ(r )i t , and h ĵðr; tÞi t , i.e., h Ôðr; tÞi t ¼ h Ôi½rð0Þ; Uðr i ;r j Þ; Kðt; t 0 Þ; hnðrÞi t ; h ĵðr; tÞi t ðr; tÞ .
An important consequence from the one-to-one character of G is that knowledge of both hnˆ(r )i t and h ĵðr; tÞi t is equivalent to the knowledge of r(t) itself (while holding r(0), U(r i ,r j ), and K(t,t 0 ) fixed).Therefore, if we accomplish to reproduce hnˆ(r )i t and h ĵðr; tÞi t via an alternative procedure other than (1) and provided we have knowledge of the G map, we are in principle capable of computing all the properties hO ˆ(r,t)i t of the system since h Ôðr; tÞi t ¼ I 0 FG À1 ðhnðr Þi t ; h ĵðr; tÞi t Þ.Given this context, we proceed to prove that G is injective.
Theorem 1. G is a one-to-one map.Proof.We generalize the original construction for closed systems by Vignale 22 and its Markovian adaptation to open systems by DADV. 16We shall present the whole proof for the sake of clarity, in spite of the fact that several steps might be repetitive to the reader who is familiar with the works cited.
We consider two systems: The original one, described by the density matrix r(t), varies in time according to (1) with the Hamiltonian H ˆS(t) given by (3); the auxiliary one, associated with the density matrix r 0 (t) and starting as r 0 (0) = r(0), evolves under the master equation in such a way that it renders the same particle and current densities as the original system at all points in time and space: We denote the expectation values taken with respect to r 0 (t) by hi 0 t .In this proof, we specialize to the case where the memory kernels in both systems are the same: K 0 (t,t 0 ) = K(t,t 0 ).
The difference between the evolution of eqn (1) and its primed analogue ( 4) is that we allow for the possibility that r(t) a r 0 (t) because although the memory kernel K(t,t 0 ) is assumed to be the same in both systems, we consider a Hamiltonian H ˆ0S (t) for the auxiliary system which may differ from its original counterpart H ˆS(t) by the vector potential: Here we also consider the particular scenario where the interparticle potential of the two systems is the same, U 0 (r i ,r j ) = U(r i ,r j ).Note that we have primed the auxiliary system's current density operator because it depends explicitly on Ã0 (r,t).We now show that in order for hypotheses (5) to be fulfilled, it must happen that Ã0 (r,t) = Ã(r,t), therefore proving the injective character of G. 35 We begin by writing the equation of motion for h ĵðr; tÞi, which can be systematically carried out by using eqn (1): Every term in this equation has a physical interpretation: q Ã(r,t)/qt and r Â Ãðr; tÞ are proportional to the electric and magnetic fields respectively, so their joint contribution represents the Lorentz force.Dðr; tÞ ¼ À is the divergence of the stress tensor, which describes the exchange of momentum between the directions a and b, where a,b = x,y,z.Fðr; tÞ is the internal force density caused by the pairwise potential Fðr; tÞ ¼ À and finally the dissipative part of the evolution is given by Gðr; tÞ ¼ Trf ĵðr; tÞð R t 0 dt 0 Kðt; t 0 Þrðt 0 ÞÞg.Similarly, for the auxiliary system, we can analogously write the equation of motion for h ĵ 0 ðr; tÞi 0 t : where the primed variables have their respective obvious meanings.By subtracting eqn ( 8) from ( 7) and imposing hypotheses (5) we can easily arrive to: where we have defined the auxiliary system vector potential with the original vector potential as the reference, Ã(r,t) = Ã(r,t) + D Ã(r,t).Therefore, D Ã(r,t) is the unknown function we are aiming to find.At this point it may seem difficult to solve for D Ã(r,t) directly from eqn (9).However, progress can be made if one assumes that the time-dependent functions are all analytic and therefore admit a Taylor expansion about t = 0.If we denote the Taylor coefficients by O k 1 k! @ k Oðr;tÞ @tk t¼0 , once we collect the terms of order t l we end up with the following identity: Note that all the Taylor coefficients are still position dependent.We now claim that the right hand side of eqn (10)  does not contain any term D Ãk (r ) for k 4 l.This is obvious for the terms that contain D Ã(r,t) explicitly.Let us study the other terms.By definition, Dl 0 ðrÞ can be written as: We know that D Ãk (r ) is the only coefficient of D Ã(r,t) that appears in @ k v0 ia @t k , whereas the highest order coefficient in @r 0 ðtÞ @t is D Ãk (r ) (see eqn (1)).Based on these two observations, it is straightforward to see that D0 l ðr Þ contains terms D ÃkÀ1 (r ) with k at the most being l.Similar arguments allow us to conclude the same k r l condition for F0 l ðr Þ.Finally, for the term G0 l ðr Þ we have: Here, the dependence of @ k ĵ 0 ðr; tÞ , we encounter two cases: In the first case, K 0 (t,t 0 ) is either a smooth function in t and t 0 or can be approximated as such, so we apply Leibniz rule k times: If we evaluate this derivative at time t = 0, the first term in the right vanishes.By examining the second term we conclude that, by analogous reasons to the ones presented above, this term only contains D ÃkÀ2 (r ) as its highest order coefficient.In summary, for the first case, every D Ãk (r ) in Gl ðr Þ has k r l À 2.
In the second case, where K 0 (t,t 0 ) is not an analytic function in t and t 0 , it is not possible to anticipate a general conclusion.However, the reader can easily check that as long as the integral R t 0 dt 0 K 0 (t,t 0 )r 0 (t 0 ) is analytic and does not contain second derivatives of r 0 (t) with respect to t or higher, the conclusion about the order of the coefficients D Ãk (r ) in Gl ðr Þ still holds.For example, in the KL form, although Lþ j Li ĵ 0 ðr; tÞ À 1 2 ĵ 0 ðr; tÞ Lþ j Li i are both analytic in t, and in fact, Gl ðr Þ only contains terms D Ãk (r ) with k r l.We believe that physical memory kernels are either smooth as in the first case, or if they belong to the second case, they follow the conditions we have outlined, in which case, we have finally proved the claim in general.
The result above allows to regard eqn (10) By hypothesis, r 0 (0) = r(0), so D Ã0 (r ) = 0; however, eqn (10)  shows that given this initial condition, D Ãk (r ) = 0 for all k, which means that Ã0 (r,t) = Ã(r,t) for all t A [0,t c1 ], where t c1 is the largest time for which the analyticity of the functions is ensured.However, the whole argument may be repeated by expanding the time dependent functions about t = t c1 until the radius of convergence in this interval is exhausted at t = t c2 , 36 and so on for the subsequent intervals, thus eventually proving that Ã0 (r,t) = Ã(r,t) for all times.This means that we have finally showed the claim of Theorem 1. Obviously, this argument fails in the pathological case where there is a vanishing radius of convergence, in which case we cannot iterate the procedure.However, just as in ref. 16 and 22, we work under the assumption that the equations of motion we are considering are smooth enough to rule this case out, although we shall acknowledge that a more rigurous study of this claim for the future would be of interest.Finally, one could imagine extending this proof to piecewise analytic vector potentials Ã(r,t).We refer the reader to ref. 37 for more details on this regard.Before we close this section, we carry out a consistency check.Just as we derived the equation of motion for the current densities, we can prove from ( 1) and (4) that the particle densities evolve as: Interestingly, eqn ( 15) and ( 16) do not exhibit the continuity form we are used to, but rather have leakage terms which are related to the memory kernels K(t,t 0 ) = K 0 (t,t 0 ).We discuss the implications about this issue in the next section.However, this apparent anomaly does not interfere with the above proof.The reason is that since we have concluded that the vector potentials in both systems must be the same, Ã0 (r,t) = Ã(r,t), from eqn ( 6) and (3) the Hamiltonians must satisfy H ˆ0S (t) = H ˆS(t), and therefore, from eqn (4) and (1), we conclude that the density matrices are equal at all times, r 0 (t) = r(t).This means that eqn ( 16) and ( 15) are identical, ensuring the consistency of our proof.Although the points we are making in this paragraph might seem trivial, in the next sections, they will be quite important for the establishment of a KS scheme for TD-CDFT.Notice that as opposed to standard TD-CDFT, the one-toone map G we have established goes from Ã(r,t) to the combination ðhnðr Þi t ; h ĵðr; tÞi t Þ as opposed to just h ĵðr; tÞi t .The reader might wonder if there is any redundance in our definition of the G map; for instance, we could imagine just defining G: Ãðr; tÞ !h ĵðr; tÞi t instead of considering the codomain to be the set of Cartesian products ðhnðr Þi t ; h ĵðr; tÞi t Þ.The map would in fact be redundant if the continuity equation were satisfied: @hnðr Þi t @t ¼ À r Á h jðr; tÞi t , which is true for closed systems.The reason is because knowing the function h ĵðr; tÞi t would imply knowing @hnðr Þi t @t , but since by hypothesis we know r(0), we can calculate hnˆ(r )i 0 and thus by integration over t, also compute hnˆ(r )i t for all r and t.However, as shown explicitly in eqn ( 15) and ( 16), continuity does not hold in general in the context we are working on, so the claims we are making in this article about the G map would not necessarily be true if we consider h ĵðr; tÞi t alone, rather than the combination ðhnðr Þi t ; h ĵðr; tÞi t Þ, as the codomain of G.

On the problem of leakage and the establishment of a Kohn-Sham scheme
The fact that ( 15) is not a continuity equation was pointed out several years ago by Frensley, 38 but has only recently been extensively studied independently by Gebauer and Car (GC) 39 and later on by Bodor and Diosi (BD). 40These authors have focused on the case where a KL master equation is derived microscopically from a system potential and a system-bath interaction both of which are local in the system coordinates.Let us briefly comment on their work in light of our present study.
First of all, let us consider the total object composed of the system S and the bath B with Hamiltonian H ˆ= H ˆS + H ˆB + H ˆSB .Here SB denotes the system-bath interaction.If H ˆSB and the potential in H ˆS depend locally on the coordinates of the system, it can be easily shown that the exact particle and current densities (obtained from the unitary evolution of the composite density matrix) do satisfy the continuity equation and the leakage term vanishes.However, depending on the level of description of the bath, which is associated with the memory kernel K(t,t 0 ), the leakage term Tr{nˆ(r,)( R t 0 dt 0 K(t,t 0 )r(t 0 ))} might not necessarily be equal to zero.In fact, for the particular case that GC and BD have studied, which is the KL master equation derived microscopically from H ˆ, the leakage term becomes where the jump operators L ˆi are in general non-local in space.Therefore, nˆ(r,t) does not necessarily commute with these operators and each of the traces in the right hand side of eqn (17) might not vanish.Physically, leakage appears in the KL formalism due to the imposition of the Born-Markov approximation, which coarse-grains the events that occur at the time-scale of the bath (see ref. 40).The object h ĵðr; tÞi t as calculated naively by taking the expectation value of the current density operator h ĵðr; tÞi with respect to the density matrix r(t) evolved through eqn (1) with the memory kernel given by ( 2 the contribution that arises from the instantaneous redistribution of position and momenta of the system due to collisions with the bath.Therefore, the total current is composed of h ĵðr; tÞi t , which is the object we have been studying throughout this article and which GC and Piccinin (GCP) 41 have called the Hamiltonian current, plus the non-Markovian feedback we just described, which has been named the collision current.It is noteworthy, however, that the authors were able to derive an explicit expression for the collision current operator ĵcol ðr; tÞ which can be used to calculate the collision current by taking the trace of its action on the (coarse-grained) density matrix r(t) given by eqn (1), that is, by computing the quantity h ĵcol ðr; tÞi.In other words, it is still fine to use eqn (1) but we need to keep in mind that if we are asked to calculate the total current density we will need to compute h ĵðr; tÞ þ ĵcol ðr; tÞi t rather than just h ĵðr; tÞi t .On the other hand, h ĵðr; tÞi t is not void of physical meaning, since it is associated to the work done by the flow of particles. 41From a more pragmatic perspective, even if we disregard the physical meaning of h ĵðr; tÞi t for a moment, Theorem 1 still guarantees us that the knowledge of hnˆ(r,t)i t and h ĵðr; tÞi t for all r and t amounts to the same information as r(t) for all t, which in turn implies knowledge of h ĵcol ðr; tÞi t , and therefore of the total real current.Obviously if instead of staying in the KL form, we use the exact memory kernel K(t,t 0 ) derived microscopically from H ˆ, the continuity equation would be fully restablished.For the intermediate case of non-Markovian memory kernels, it is not very clear at this moment how would the ''leakage'' in eqn (15)  be affected and what h ĵðr; tÞi t would precisely mean, though these are definitely interesting questions to explore.Going beyond the case studied by GCP, in many situations of practical interest, master equations used to model relaxation from non-equilibrium configurations will contain memory kernels K(t,t 0 ) that do not necessarily cause the rightmost term in eqn (15) to vanish, in which case, an appropriate interpretation for the object h ĵðr; tÞi t is desirable.From a more general perspective, the claim that Theorem 1 provides suffices to regard h ĵðr; tÞi t as a quantity that has computational merit on its own, and this pragmatic approach will turn out to be good enough to work through the rest of this article.Hereafter, unless otherwise specified, the word current will be used again to refer to h ĵðr; tÞi t since we will not worry about the other types of currents.
After this interlude, the reader might question why should we care about leakage if our version of the RG theorem did not seem to be affected by it.However, let us ask the following question: Can we suggest a theorem whereby the same particle and current densities hnˆ(r )i t and h ĵðr; tÞi t of the original system can be reproduced by a non-interacting open KS system with the same memory kernel K(t,t 0 ) and starting with a density matrix r 0 (0) which may be different from r(0)? DADV have argued that this is possible for the case of Markovian master equations. 16If we follow their rationale, we just need to make the following trivial modifications to the proof of Theorem 1: (a) Set U 0 (r i ,r j ) = 0 in (6), and (b) lift the restriction that r 0 (0) = r(0).At a first glance, it may seem that such a proof is valid.Nevertheless, the problem manifests itself when we arrive to the consistency check in (15) and (16).Let us analyse the situation carefully: In constructing the KS vector potential Ã0 (r,t), we considered both the hypotheses in eqn ( 5) and the set of eqn ( 7) and ( 8), which describe the motion for h ĵðr; tÞi t and h ĵ 0 ðr; tÞi 0 t respectively.However, the constraints given by eqn ( 15) and ( 16), were not used at all, but should be satisfied by the respective particle and current densities.The only way the constructed Ã0 (r,t) does what we want (which is to reproduce hnˆ(r,t) t and h ĵðr; tÞi t from the original system in the auxiliary system, as expressed in ( 5)) is if the ''leakage'' terms in both systems are the same: If K 0 (t,t 0 ) = K(t,t 0 ), ( 18) becomes an additional restriction which cannot be satisfied in general.DADV have suggested that these terms normally vanish; 42 however, as we have commented before, based on the work of GCP 39 and BD, 40 these terms are in general not negligible: In fact, if the system-bath interaction is local in the coordinates of the system, from the microscopic derivation of the KL master equation, it can be seen that L ˆa is a non-local function in space, which then does not commute with nˆ(r ), and hence, Tr{nˆ(r )L(r(t)) does not necessarily vanish).Interestingly however, due to the weaker conditions expected for TD-DFT, the analogous KS scheme suggested by BCG 14 is actually justified.For a discussion of this case, we refer the reader to ref. 17.
The reason why Vignale did not encounter this problem for a closed system (K(t,t 0 ) = 0) in (ref.22) is because in that case, there are no ''leakage'' terms and the equations ( 15) and ( 16) are just redundant with respect to (5).It seems, however, that the only problem in DADV's scheme is a matter of the restriction of two variables to reproduce, but only constructing a single function.We might imagine that lifting some of the assumptions in the suggested KS scheme will give freedom to the mathematical structure and will solve the problem.In fact, in the next section, we show that by not assuming that K 0 (t,t 0 ) must equal to K(t,t 0 ), a KS scheme is actually possible.Additionally, in section 5, we will prove that there is a simpler alternative KS scheme where the auxiliary system is taken to be both closed and non-interacting.

A possible modification of D'Agosta and Di Ventra's Kohn-Sham scheme
Without more preamble, using the same language that we have been employing before, we make the following claim: Theorem 2. It is possible to reproduce both the particle and current densities of the original many-body interacting OQS (hnˆ(r )i 0 t = hnˆ(r )i t , h ĵ 0 ðr; tÞi 0 t ¼ h ĵðr; tÞi 0 t for all r and t) with an auxiliary Markovian OQS with interparticle potential U 0 (r i ,r j ) starting in the state r 0 (0) at the expense of a given vector potential Ã0 (r,t) and a memory kernel K 0 (t,t 0 ) of the KL form with a single jump operator: Comments.(a) As opposed to the proof for Theorem 1, here we do not assume that in general r 0 (0 non-interacting auxiliary system, i.e.U 0 (r i ,r j ) = 0 corresponds to the KS case.(c) Note that the theorem claims that even if the kernel of the original system K(t,t 0 ) is non-Markovian, the kernel in the KS system can be taken to be Markovian.(d) For the case where the original system's kernel K(t,t 0 ) is of the KL form, Theorem 1 becomes a fix for DADV's KS scheme.Proof.Compared to the proof given for Theorem 1, we need to construct both Ã0 (r,t) and K 0 (t,t 0 ) instead of just Ã0 (r,t).We follow the same steps as in Theorem 1, except that we consider the adaptations indicated in Comment (a).When we reach eqn (10), we still follow all the arguments about how there are no terms D Ãk (r ) for k 4 l in the right-hand side.However, we do not proceed further to solve for each of the coefficients D Ãl+1 (r ) since G0 l ðr Þ contains K 0 (t,t 0 ), which at the moment is an unknown function since we have not constructed it.Having this in mind, we combine eqn ( 10) and ( 18) in order to solve for both D Ã(r,t) and K 0 (t,t 0 ).There could be many possibilities for K 0 (t,t 0 ), but as we have specified in the claim, we choose the simplest, namely, the KL form K 0 (t,t 0 )r(t 0 ) = L 0 (r(t 0 ))d(t À t 0 ) with the respective generator L ˆ0+ (t).
Let the leakage term be denoted by l(r,t) Tr{nˆ(r,t) R t 0 K(t,t 0 )r(t)} so that ( 18) can be read as l(r,t) = l 0 (r,t).The l-th time derivative of this equation evaluated at t = 0 reads l l (r ) = l 0 l (r ).Here, we employ the same notation for the Taylor expansion coefficients as before.This equation yields a possible solution (there are many) for a matrix L 0 l in terms of the appropriate time derivatives of Ã(r,t), K(t,t 0 ), ĵðr; tÞ, r(t), r 0 (t), and lower order derivatives of L ˆ0 at t = 0.By similar arguments as the ones presented in section II, one can show that the chosen solution for L 0 l will only contain terms D Ãk (r ) such that k r l.On the other hand, the recursion of eqn (10) shows that the first step in which L 0 l appears is on the right hand side of the solution for D Ãl+1 (r ) (inside of G0 l ðr Þ, see eqn (10)).With these observations in mind, the construction of the coefficients L 0 l and D Ãl (r ) for all l Z 0 is possible through a small modification of the recursion (10).By using strong induction, we assume that right before the lth step of the recursion, we already know D Ãk (r ) and L k for all k o l.The lth step consists of (a) the solution of l l (r ) = l 0 l (r ) for L l , followed by (b) the solution of (10) for D Ãl+1 (r ).Symbolically: Á Á Á-D Ãl (r ) -L ˆl -D Ãl+1 (r ) -Á Á Á.Note that D Ã0 (r ) is still given by ( 14), and in fact, once we know its value, we can solve l(r,0) = l 0 (r,0) (which depends on the already obtained D Ã0 (r )) for a possible L 0 , therefore building a base for the recursion.This concludes the inductive proof.
With Theorem 2, we have not only provided a fix for DADV's KS scheme for Markovian OQS, but we have also generalized it for their non-Markovian counterpart.Although this scheme might be advantageous in some practical applications, we believe that the alternative scheme we propose in the next section offers a new perspective in the formulation of KS schemes and might be more practical than the one we considered here.

An alternative KS scheme for open systems
The next theorem constitutes our alternative proposal for KS scheme to be used in the practice of TD-CDFT for an arbitrary OQS: Theorem 3. (A,C Kohn Sham scheme) There exists an auxiliary closed KS system which starting in the state r KS (0) evolves unitarily under a KS Hamiltonian H ˆKS (t): producing particle and filtered current densities that are related to the particle and current densities of the original many body interacting open quantum system governed by (1) by: where we have defined the filtered current operator as ĵfilt ðr; tÞ ¼ P i 1 2 f vi;filt ð ri ; tÞ; dðr À ri Þg with the filtered velocity operator given by vi;filt ðr i ; tÞ ¼ pi þe ÃKS ð ri ;tÞ m .Additionally, we denote the expectation values in the KS system in the usual form hO ˆ(r,t)i KS,t = Tr(O ˆ(r,t)r KS (t)).We shall call ÃKS (r,t) the KS vector potential, and C(r,t) the leakage potential.
Proof.We proceed similarly to the proof in Theorem 1: Given the master eqn (1), we evolve it to compute hnˆ(r )i t and h ĵðr; tÞi t for all positions and times.We explicitly ''solve'' for ÃKS (r,t) and C(r,t) such that the density matrix where T denotes the time-ordering operator, renders the expectation values hnˆ(r )i KS,t and h ĵfilt ðr; tÞi KS;t which satisfy the relations given in eqn (21).The explicit construction of ÃKS (r,t) and C(r,t) implt their existence, and therefore prove the theorem.
We construct these potentials in two steps.
Step 1. (Construction of C(r,t)).The equation of motion for the particle density of the open system is just the pseudocontinuity eqn (15).On the other hand, for the KS system, we can write: By applying the constraints (21), it follows that C(r,t) must be given by where the addition of any function in time as an integration constant is valid at this point.It is possible that a particular way of choosing boundary conditions may enhance the intuitive physical meaning of the leakage potential C(r,t).
For instance, for the particular case studied by Gebauer, Car, and Piccinin, we can choose C(r,t) so that e Cð r;tÞ m hnðr Þi KS;t is equal to what they have termed the collision current 41 (see discussion in Section III).However, this choice is not formally necessary and C(r,t) can simply be regarded as a strategy that reestablishes the continuity equation via Step 2. (Construction of ÃKS (r,t)).The equation of motion for h ĵðr; tÞi t for the original system was given in (7).For the closed KS system, we can analogously write an equation of motion for h ĵðr; tÞi KS;t : @h ĵðr; tÞi KS;t @t ¼ hnðr Þi KS;t m @ð ÃKS ðr; tÞ þ Cðr; tÞÞ @t Note the absence of a dissipative GKS ðr; tÞ term due to the unitarity of the evolution of the KS system.Combining the constraints ( 21) with ( 7) and ( 24) yields: where as usual we have defined Ã = ÃKS (r,t) + D Ã(r,t).Just as in (10), we write out the terms of (25) of order t l : From analogous arguments to the ones in the proof for Theorem 1, all the terms D Ãk in the right hand side of ( 26) are such that k o l.This is a consequence of the linearity of the master eqn (1).The recursion ladder can be solved if we know the value of D Ã0 (r ), which from the constraints (21)   for t = 0 yield D Ã0 ðr Þ ¼ Note that provided we decide on a systematic way to choose the boundary conditions for C(r,t), Ã(r,t) is uniquely defined.QED.
We call the reader's attention to several features of Theorem 2: INTERACTING TO NON-INTERACTING.Just as in the regular KS schemes, we can set the auxiliary system to be non-interacting, U 0 (r i ,r j ) = 0.
MIXED STATES TO PURE STATES.For practical applications, it might be the case that sometimes it is easier to work with pure states rather than mixed states.In fact, from a computational perspective, the numerical propagation of a pure state in an R-dimensional Hilbert space requires the propagation of 2(R À 1) real parameters of a vector, whereas for a mixed state, it requires the evolution of R 2 À 1 real parameters.If the KS initial state is taken to be pure, i.e. r KS (0) = |c KS (0)ihc KS (0)|, where |c KS (0)i is the initial KS wavefunction, Theorem 2 guarantees that the KS state will stay pure at all times due to the unitary evolution via There is however, the subtle question of whether given any master eqn (1), |c KS (0)i would always exist such that it can satisfy the constraints (21) at t = 0.However, the first constraint with respect to the particle densities can be trivially satisfied for bosons, and by the so-called Harriman construction for the case of fermions, which is a prescription to construct an N particle single Slater determinant with a specified particle density. 44nce this initial wavefunction |c KS (0)i is constructed, the constraint in the current density can always be satisfied by choosing D Ã(r,0) as explained at the end of the proof for Theorem 2.
OPEN SYSTEM TO CLOSED SYSTEM.Theorem 2 provides us with the formal footing to calculate the particle and current densities of an open system in a closed system.The artifact involved in this mapping is the construction of the leakage potential C(r,t) and the definition of a filtered current operator ĵfilt ðr; tÞ.If there is no leakage, C(r,t) vanishes identically for all r and t.This happens naturally for K(t,t 0 ) = 0, which is the case of closed systems, but might also happen in open systems in the special case where the system-bath interaction is local in the coordinates of the system.In such case, ĵfilt ðr; tÞ ¼ ĵðr; tÞ and we effectively have a map from an open to a closed system without any additional artifact.However, as discussed in section 3, this case might be an exception rather than a rule.In any case, the use of a closed system as the auxiliary KS system might turn out to have very important implications, since at the present time, there are many more numerical propagation schemes for closed quantum systems, 43 than their open counterpart.
An important feature associated with the use of a closed KS system is that we do not need to worry about the propagation of the Slater determinant as a whole, but rather evolve each orbital unitarily and reconstruct the Slater determinant at all points in time when needed, just as one might envision performing real time TD-CDFT of closed systems.This observed numerical consistency with the inverted wavefunction c KS (r,t).Going back to the gauge where the scalar potential vanishes, by recognizing that the total vector potential ÃKS (x,t) + C(x,t) in eqn (20) must equal to x @ @x R t 0 V KS ðx; t 0 Þ dt 0 , we finally constructed ÃKS (x,t) by performing another step of trapezoidal rule numerical integration and substracting out the already constructed leakage potential C(x,t).To gain more intuition on the problem, we have plotted c KS (r,t), rather than the gauge-transformed wavefunction cKS (r,t) = e ikx 2 t/2 c KS (r ) since the latter exhibits slightly more complicated beating patterns than the more familiar form shown by c KS (r,t) which was reconstructed in the gauge where the external potential is purely scalar.
We discuss the results by referring to the plotted quantities.In order to have a better visualization for the KS scalar potential V KS (r,t), we plotted the function with a time dependent shift s(t): V KS (r,t) -V KS (r,t) + s(t), so that V KS (0,t) + s(t) = 0 for all t, which obviously does not change the dynamics of the system.For the case of no dissipation (z = 0) we encounter the expected results (Fig. 2): Rabi oscillations in time for the particle (2a) and current densities (2b) as well as for the wavefunction (2c,2d).The current density does not vanish because at every time point, the state of the system is a superposition of the eigenstates |3i and |4i.The reconstructed KS scalar potential V KS (r,t) fortunately yields the expected harmonic well for all times (2e).Since the KS vector potential ÃKS (x,t) is reconstructed from V KS (r,t), it exhibits the expected form ÃHO (x,t) = xˆkxt (2f).The result for the leakage potential C(r,t) is theoretically zero for all r and t but in our calculations, there is a noise on the order of B10 À4 , which compared to the energy scale of order B10-100, can be neglected (2g).
For both cases with dissipation (Fig. 3 and 4), the Rabi oscillations of the particle and current densities are smoothed in time (3a, 3b, 4a, and 4b).At long times, hnˆ(x)i t becomes the Gaussian function associated with the ground state |0i, and h ĵðx; tÞi t vanishes for all x.The KS wavefunction for the z = 0.15 case resembles the one for the situation z = 0 although its regular oscillatory pattern is slowly smoothed to become the ground state of the harmonic oscillator.Therefore, the real part slowly becomes a Gaussian function as time increases (3c), whereas the imaginary part simply decays (3d).For z = 0.5, the example with strong dissipation, we see that the same effects for the KS wavefunction but with higher intensity.The real part of the KS wavefunction quickly becomes Gaussian (4c) just as the particle density (4a), and the corresponding imaginary part vanishes only after a short period of time.The leakage potentials C(r,t) in both cases (3g and 4g) are relatively small and are only present during the first periods of the evolution of the system.This observation might not be true for other types of dissipation.It is interesting to note that V KS (r,t) for z = 0.15 (3e) differs significantly from its non-dissipative counterpart (2e).It can be qualitatively described as a series of multi-well potentials that vary in position as time progresses in order to accomodate most of the oscillatory particle density.However, for the z = 0.5 case, V KS (r,t) oscillates strongly only for early times but afterwards, it becomes $ mo 2 x 2 2 .The latter observation cannot be perceived very clearly from Fig. 4e due to the length scale of the z-axis.However, the fact that the KS vector potential slowly becomes parabolic is consistent with the fact that the particle density simultaneously becomes Gaussian.Finally, whereas the constructed vector potential ÃKS (r,t) for the z = 0.15 case (3f) differs significantly from the z = 0 case (2f), the analogous plot for the z = 0.5 case (4f) exhibits a structure that reminds the planar form of its non-dissipative counterpart.This is a consequence of the similarity between the scalar potentials of the z = 0 and z = 0.5 cases, and the null and negligible intensities of the leakage potentials for these scenarios.

Conclusions
In this article we have systematically developed a TD-CDFT formalism for generalized many-body open quantum systems.We extend the previous work on this field (ref.14 and 16) to the non-Markovian regime, address issues of representability which have not been discussed before, and suggest novel ways to conceive KS schemes for non-Hamiltonian dynamics.In the next few paragraphs, we enumerate the main results of our investigation, comment on their relationship with previous studies in the field, and finally also indicate some open questions that remain to be addressed in our future work.
In section 2, we proved the RG theorem for TD-CDFT in the context of generalized OQS.To do so, we adapted the construction given by Vignale in ref. 22 for closed systems.For the special case of Markovian dissipation of the KL form, our result reduces to one by DADV in ref. 16.We can highlight two main differences between our version of the RG theorem and the one given by Vignale for TD-CDFT for closed systems.First, our statement in Theorem 2 requires an additional specification for the G map, namely fixed memory kernel K(t,t 0 ), besides fixed initial density matrix r(0) and interparticle potential U(r i ,r j ).Second, we have defined the image set of G to consist of coordinate pairs of particle and current densities ðhnðr Þi t ; h ĵðr; tÞi t Þ instead of just particle densities h ĵðr; tÞi t , as one would have done in the original formulation of TD-CDFT.It turns out that this specification is not redundant but necessary, since specifying h ĵðr; tÞi t does not specify hnˆ(r )i t due to the lack of continuity in master equations (eqn ( 15) and ( 16)).
In section 3, we analysed the currently proposed KS scheme for TD-CDFT for OQS given by DADV and we concluded that there is a slight inconsistency that renders the scheme non-applicable to many cases.However, in section 4 we have suggested a possible way to restore the method, where we do not only need to construct a KS vector potential, but also a KS dissipation, which we have claimed can be taken to be Markovian in general.However, with the goal of envisioning an alternative to fix this problem, in section 5 we suggested a novel KS scheme where the auxiliary system is both non-interacting and closed.Theorem 3, which we call the A,C KS scheme, constitutes the main result of our article.In a practical level, such as with the aim of carrying out numerical simulations, this property might turn out to be very practical since the quantum mechanical propagation of a closed system is numerically simpler and more efficient than the analogous exercise for OQS.Also, by using a closed KS system, we avoid a computational nuisance which is present in the schemes of BCG and DADV, namely, that the real-time propagation of the Slater determinant via a master equation or a stochastic Schro¨dinger equation, requires non-unitary jumps between different orbitals, which does not allow for the independent evolution of each orbital.Our proposal involves the independent unitary evolution of each orbital, therefore resembling the real-time propagation of orbitals in regular TD-CDFT.
Our future work on this subject will focus on the implementation of functional approximations for ÃKS (r,t) and C(r,t).Just as with the Vignale-Kohn functional, 24 we might regard the electron liquid as our starting model system.However, in our case we must consider the effects that dissipative dynamics have on the exchange and correlation properties of the particles in the system.In section 4, we have enumerated a series of properties that these functionals need to satisfy; we shall use these properties as guidelines for the development of approximate functionals.Ultimately, our goal is to apply this formulation to a realistic system where the non-Hamiltonian evolution is essential to understand its function, such as excitation transfer in chromophores embedded in a protein environment 10 or resonant surface enhanced Raman scattering. 4 In summary, our work provides an alternative tool to study generalized open quantum systems within the realm of TD-CDFT.

Fig. 1
Fig.1Scheme showing the different maps that are relevant to the proof of the Runge-Gross (RG) theorem for generalized open quantum systems.By holding the initial density matrix r(0), the interparticle potential U(r i ,r j ), and the memory kernel of the dissipation K(t,t 0 ) fixed, we define the map G goes from vector potentials Ã(r,t) to particle and current densities ðhnðrÞi t ; h ĵðr; tÞi t Þ.Our adaptation of the RG theorem for open quantum systems claims that G is one-to-one, which implies that it is possible to define an inverse map G À1 from the image set of G (see sections 1 and 2).Here, ME stands for master equation.
as a recursion relation to compute all the coefficients D Ãk (r ) provided we know D Ã0 (r ).To obtain D Ã0 (r ), we recognize that This journal is c the Owner Societies 2009 Downloaded by Harvard University on 02 March 2012 Published on 11 May 2009 on http://pubs.rsc.org| doi:10.1039/B903064FD Ã0 (r ) = D Ã(r,0) and that the constraints of eqn (5) at time t = 0 yield: ) might not necessarily correspond to the total measured current in a nanodevice, since it misses This journal is c the Owner Societies 2009 Phys.Chem.Chem.Phys., 2009, 11, 4509-4522 | 4513 Downloaded by Harvard University on 02 March 2012 Published on 11 May 2009 on http://pubs.rsc.org| doi:10.1039/B903064F r j ).(b) The corollary of Theorem 2 for the 4514 | Phys.Chem.Chem.Phys., 2009, 11, 4509-4522 This journal is c the Owner Societies 2009 Downloaded by Harvard University on 02 March 2012 Published on 11 May 2009 on http://pubs.rsc.org| doi:10.1039/B903064F

Fig. 3 1
Fig.31-D harmonic oscillator coupled to a heat bath that induces decays to the ground state with a rate z = 0.15.Just as in Fig.1, the particle (a) and current densities (b) were computed through the evolution of a density matrix starting in the pure state r(0) = 1 2 (|3i + |4i)(h3| + h4|).The real part of the constructed Kohn-Sham wavefunction (c) slowly approaches the Gaussian function form as time progresses.The corresponding imaginary part (d) gradually vanishes as expected.The inverted scalar potential (e) can be qualitatively described as a series of double well potentials with time varying minima which accommodate most of the particle density in them (see (a)).The constructed vector potential (f) differs considerably from its counterpart in Fig.1.The intensity of the leakage potential (g) is very weak for this system and type of dissipation and has non-negligible values only at the beginning of the evolution.

Fig. 4 1
Fig.41-D harmonic oscillator coupled to a heat bath that induces decays to the ground state with a rate z = 0.5.Just as in Fig.1 and 2, the particle (a) and current densities (b) were computed through the evolution of a density matrix starting in the pure state r(0) = 1 2 (|3i + |4i)(h3| + h4|).The real part of the constructed Kohn-Sham wavefunction (c) quickly approaches the Gaussian function form due to the strong dissipative effect of the bath.The corresponding imaginary part (d) quickly vanishes as expected.The inverted scalar potential (e) oscillates strongly at short times, but quickly settles in the expected parabolic form.Because of the latter, the constructed vector potential (f) shares a similar tilt as its counterpart in Fig.1.The intensity of the leakage potential (g) is also weak as in Fig.(2).

1 2
Pi fdðr À ri Þ; vi ð ri ; tÞg where the canonical velocity operator vi ð ri ; tÞ is explicitly time dependent via the vector potential: vi ð r; tÞ ¼ 1 m ð pi þ e Ãð ri ; tÞÞ.The expectation value for an arbitrary observable of the system O ˆ(r,t) will be computed as usual by taking the trace with the density matrix r(t), hO ˆ(r,t)i t = Tr(O ˆ(r,t)r(t)), where hi t This journal is c the Owner Societies 2009 Downloaded by Harvard University on 02 March 2012 Published on 11 May 2009 on http://pubs.rsc.org| doi:10.1039/B903064Fhnˆ(r )i KS,t .In fact, notice that the actual (unfiltered) KS current density operator ĵðr; tÞ ¼ P Þg where vi ðr; tÞ ¼ pi þeð ÃKS ð ri ;tÞþ Cð ri ;tÞÞ m does satisfy the continuity equation by construction since ĵðr; tÞ ¼ ĵfilt þ e m Cðr; tÞnðr Þ.