Open Quantum Systems: Density Matrix Formalism and Applications

In its original formulation, TDDFT addresses the isolated dynamics of electronic systems evolving unitarily (Runge and Gross 1984). However, there exist many situations in which the electronic degrees of freedom are not isolated, but must be treated as a subsystem imbedded in a much larger thermal bath.


Introduction
In its original formulation, TDDFT addresses the isolated dynamics of electronic systems evolving unitarily [Runge 1984].However, there exist many situations in which the electronic degrees of freedom are not isolated, but must be treated as a subsystem imbedded in a much larger thermal bath.The theory of open quantum systems (OQS) deals with precisely this situation, in which the bath exchanges energy and momentum with the system, but particle number is typically conserved.Several important examples include vibrational relaxation of molecules in liquids and solid impurities, coupling to a photon bath in cavity quantum electrodynamics, photo-absorption of chromophores in a protein environment, electron-phonon coupling in singlemolecule transport and exciton and energy transfer nanomaterials.Even with simple system-bath models, describing the reduced dynamics of many correlated electrons is computationally intractable.Therefore, applying TDDFT to OQS (OQS-TDDFT) offers a practical approach to the many-body opensystems problem.
The formal development of the theory of OQS begins with the full unitary dynamics of the coupled system and bath, described by the Von Neumann equation for the density operator dρ(t) dt = −i[ Ĥ(t), ρ(t)] . (9.1) Here, Ĥ(t) = ĤS (t) + α ĤSB + ĤB (9.2) is the full Hamiltonian for the coupled system and bath and is the Hamiltonian of the electronic system of interest in an external potential v ext (r, t).This potential generally consists of a static external potential due to the nuclei and an external driving field coupled to the system such as a laser field.For an interacting electronic system, v ee (r, r ′ ) = 1/|r − r ′ | is the two-body Coulomb repulsion.The system-bath coupling, ĤSB , acts in the combined Hilbert space of the system and bath and so it couples the two subsystems.Typically, for a single dissipation channel, the system-bath coupling is taken to have a bilinear form, ĤSB = − Ŝ ⊗ B , (9.4) where B is an operator in the bath Hilbert space which generally couples to a local one-body operator Ŝ = N i=1 ŝ( pi , ri ) in the system Hilbert space.Implicit in OQS is a weak interaction between the system and bath, so that one can treat the system-bath coupling perturbatively by introducing the small parameter α as in Eq. (9.2).ĤB is the Hamiltonian of the bath, which will typically consist of a dense set of bosonic modes such as photons or phonons.The density of states of ĤB determines the structure of reservoir correlation functions, whose time-scale in turn determines the reduced system dynamics.
The goal of the theory of open quantum systems is to arrive at a reduced description of the dynamics of the electronic system alone, by integrating out the bosonic modes of the bath.In this way, one arrives at the quantum master equation, which describes the non-unitary evolution of the reduced system in the presence of its environment.In the next section, we derive the manyelectron quantum master equation and discuss common approximations used to treat the system-bath interactions.We then formulate the master equation approach to TDDFT rigorously, by establishing a van Leeuwen construction for OQS.Next, we turn to a practical Kohn-Sham (KS) scheme for dissipative real-time dynamics and finally discuss the linear response version of OQS-TDDFT, giving access to environmentally broadened absorption spectra.

The generalized quantum master equation
We begin this section by deriving the formally exact many-electron quantum master equation using the Nakagima-Zwanzig projection operator formalism [Nakajima 1958, Zwanzig 1960, Zwanzig 2001].Using projection operators, the master equation can be systematically derived from first principles starting from the microscopic Hamiltonian in Eq. (9.2) (i.e.without phenomenological parameters).This is particularly amenable to TDDFT, in which the electronic degrees of freedom are treated using first principles as well.We will then discuss the Born-Markov approximation and the widely used Lindblad master equation [Lindblad 1976].9.2.1 Derivation of the quantum master equation using the Nakagima-Zwanzig projection operator formalism Our starting point is Eq.(9.1) for the evolution of the full density operator of the coupled system and bath, where L(t) is the Liouvillian superoperator for the full evolution defined by Eq. (9.5).It may be separated into a sum of Liouvillian superoperators as (9.6)where each term acts as a commutator on the density matrix with its respective part of the Hamiltonian.Our goal is to derive an equation of motion for the reduced density operator of the electronic system, defined by tracing the full density operator over the bath degrees of freedom.
To achieve this formally, we introduce the projection superoperators P and Q.The operator P is defined by projecting the full density operator onto a product of the system density operator with the equilibrium density operator of the bath, P ρ(t) = ρeq B ρS (t) .
(9.8) Q = 1 − P projects on the complement space.In this sense, P projects onto the degrees of freedom of the electronic system we are interested in, while Q projects onto irrelevant degrees of freedom describing the bath dynamics.
Using these projection operators, Eq. (9.5) can be written formally as two coupled equations: If Eq. (9.9b) is integrated and substituted into Eq.(9.9a), one obtains By performing a partial trace of both sides of Eq. (9.10) over the bath degrees of freedom, one arrives at the formally exact quantum master equation is the memory kernel.
arises from initial correlations between the system and its environment [Meier 1999].Equation (9.11) is still formally exact, as ρS (t) yields the exact expectation value of any operator depending on the electronic degrees of freedom.
In practice, approximations to the memory kernel and initial correlation term are needed.

The Markov approximation
One often invokes the Markov approximation, in which the memory kernel is local in time and the initial correlations vanish, i.e.
The Markov approximation is valid when τ S ≫ τ B is satisfied, where τ S is the time-scale for the system to relax to thermal equilibrium and τ B is the longest correlation time of the bath [Breuer 2002].Roughly speaking, the memory of the bath is neglected because the bath decorrelates from itself before the system has had a chance to evolve appreciably [Van Kampen 1992].The timescale τ S is inversely related to the magnitude of the system-bath coupling, and so a weak interaction between the electrons and the environment is implicit in this condition as well.The Lindblad form of the Markovian master equation, (9.16) is constructed to guarantee positivity of the density matrix.This is desirable, since the populations of any physically sensible density matrix should remain positive during the evolution.In the Lindblad equation, in addition to the Markov approximation, one performs second-order perturbation theory in the system-bath interaction [Eq.(9.4)].These two approximations in tandem are referred to collectively as the Born-Markov approximation.So far our discussion has focused on a system of interacting electrons coupled to a bath.We now turn to the formulation of OQS-TDDFT.

Rigorous foundations of OQS-TDDFT
In order to formally establish an OQS-TDDFT starting from the many-body quantum master equation in Eq. ( 9.11), we must first establish the opensystems version of the van Leeuwen construction [van Leeuwen 1999].This proves a one-to-one mapping between densities and potentials for non-unitary dynamics, as well as the existence of several different KS schemes [Yuen-Zhou 2010, Yuen-Zhou 2009].

The OQS-TDDFT van Leeuwen construction
Our starting point is the master equation of Eq. (9.11), which evolves under the many-electron Hamiltonian in Eq. ( 9.3).We may now state a theorem concerning the construction of an auxiliary system.
Theorem: Let the original system be described by the density matrix ρS (t), which starting as ρS (0) evolves according to Eq. (9.11).Consider an auxiliary system associated with the density matrix ρ′ S (t) and initial state ρ′ S (0), which is governed by the master equation and with the functional forms of K′ (t, t ′ ) and Ξ ′ (t) given.Here, is the Hamiltonian of an auxiliary system with a different two-particle interaction v ′ ee (r, r ′ ) which is also given.Under mild conditions, there exists an external potential v ′ ext (r, t) which drives the system in such a way that the particle densities in the original and the auxiliary systems are the same, i.e. n(r) ′ = n(r) .This implies that Tr{ ρS (t)n(r)} = Tr{ ρ′ S (t)n(r)} is satisfied for all times.
Proof : The method we use closely parallels the van Leeuwen construction given for unitary evolution [van Leeuwen 1999].By using Eq. ( 9.11) we can find an equation of motion for the second derivative of the particle density of the original system.This is done by first deriving the equation of motion for the particle density as well as for the current density, We then differentiate both sides of Eq. ( 9.19) with respect to time and use Eq.(9.20) to eliminate the current.One thus arrives at, i.e. we demand that the densities and their first derivatives be the same at the initial time.In Eq. (9.21), each term has a clear physical interpretation.
The quantity ∇v ext (r, t) is proportional to the external electric field acting on the system, is the divergence of the stress tensor, where α, β = x, y, z label Cartesian indices, and is the internal force density caused by the pairwise potential.In addition to these quantities which arise in usual TDDFT, we have defined two new quantities, which are unique to OQS-TDDFT and arise from forces induced by the bath.
We can now repeat the same procedure in the primed system, to arrive at the equation of motion for the second derivative of the density in terms of primed quantities, If we subtract Eq. (9.21) from Eq. (9.26) and assume that n(r) ′ t = n(r) t , we arrive at the equation where we have defined ∆v ext (r, t) ≡ v ′ ext (r, t) − v ext (r, t).We now expand both sides of Eq. (9.27) in a Taylor series with respect to time to arrive at The left-hand side of Eq. (9.28) contains Taylor coefficients of v ′ ext (r, t) of order l, while the right-hand side depends only on Taylor coefficients of v ′ ext (r, t) of order k < l and known quantities.Equation (9.28) can therefore be regarded as a unique recursion relation for constructing the Taylor coefficients of the auxiliary potential v ′ ext (r, t), once a suitable boundary condition is specified.We assume that v ′ ext l (r) → 0 sufficiently quickly as |r| → ∞ for all l.
Several different KS schemes are now evident.If one sets v ′ ee (r, r ′ ) = 0, but keeps the system open by setting K′ (t) = KKS (t) and Ξ ′ (t) = Ξ KS (t), the auxiliary system is a non-interacting, but open KS system.This is similar to the construction used in [Burke 2005c], but encompasses the non-Markovian case as well.However, one may also choose v ′ ee (r, r ′ ) = 0 and K′ (t) = Ξ ′ (t) = 0, whereby the density of the original open system is reproduced with a closed (unitarily evolving) and non-interacting KS system.

The double adiabatic connection
The content of our proof is conveniently summarized by parametrizing the auxiliary system's master equation with two coupling constants λ and β as, where Here, λ scales the electron-electron interaction and lies in the range 0 λ 1.When λ = 1, we have a fully interacting system, while when λ = 0 we have a system of non-interacting electrons.The memory kernel and initial correlations are functions of λ as well.Similarly, β scales the non-unitary terms in the master equation and lies in the range 0 β 1.When β = 1, we have a fully open system while when β = 0 the system evolves unitarily.The theorem of Sect.9.3.1 guarantees the existence and uniqueness of a potential v ′ ext (λ, β, r, t) for all λ and β, which drives the auxiliary system in such a way that the true density is obtained independent of the values of λ and β.This can be viewed as a two-dimensional extension of the usual electronelectron adiabatic connection in closed-systems TDDFT [Görling 1997].It is depicted graphically in Fig. 9.1.At the coordinate (1,1), we have the original fully interacting and open system, while at the coordinate (0,0) we have the non-interacting and closed KS scheme.Defining K′ (λ = 0) ≡ KKS (t) and Ξ ′ (λ = 0) ≡ Ξ KS as the memory kernel and initial correlations of an open system of non-interacting electrons, we see that the point (1,0) describes the open KS scheme.In the remainder of the chapter, we will focus on the points (0,0) and (1,0).However, our proof shows that any coordinate lying within the double adiabatic connection square represents a viable KS scheme.As in DFT and TDDFT, the double adiabatic connection provides a powerful tool for deriving exact conditions on OQS-TDDFT functionals and is currently being explored in more detail.

Simulating real-time dissipative dynamics with a unitarily evolving Kohn-Sham system
In the previous section, we saw that it is possible to take the KS system to be a non-interacting system evolving unitarily under a time-dependent driving field that will reproduce the particle density of the original interacting OQS.In this scheme, the KS scalar potential v ′ ext can be expressed as a sum of several contributions: Here, v ext is the same external potential as the one acting on the real system (electron-nuclei and possibly an external perturbation, such as a laser field).The electron-electron interaction is replaced by the sum of a Hartree term and a standard approximation to the exchange-correlation (xc) term v xc , such as an adiabatic functional [Furche 2002a].Finally, v bath is a new term which represents a driving field that mimics the interactions of the system with the bath.This KS scheme places electron-electron and system-bath interactions on the same footing, so real-time TDDFT computer codes could in principle be easily modified to include the dissipative effects of an environment [Castro 2006].Furthermore, by propagating orbitals instead of a density matrix as in an open Kohn-Sham scheme, the scaling of the computational time with system size is the same as for a standard closed-system propagation.
It is therefore our goal to suggest practical, yet reasonably accurate approximations to v bath .Before proceeding along this direction, our goal deserves two comments: First, it goes in line with the spirit of standard real-time TDDFT, where the many-body Hamiltonian of the original system is replaced by a KS Hamiltonian which, although non-interacting, is a functional of the particle density, and therefore is propagated as a nonlinear Schrodinger equation.Second, equations of motion for classical systems coupled to heat baths are often described by Langevin equations, where frictional forces are functions of the velocity of the particle itself (Stokes law) [Zwanzig 2001], making their numerical propagation nonlinearly dependent on system variables such as the current and density.Therefore, the search for approximations to v bath could start by investigating work already done in the field of dissipative non-linear Schrödinger equations (NLSE) [Kostin 1972, Kostin 1975, Bolivar 1998, Haas 2010], as well as of time-dependent self consistent field (TD-SCF) methods [Makri 1987, López 2010, Martinazzo 2006].In this section, we describe a simple Markovian bath functional inspired by the NLSE suggested by Kostin [Kostin 1972].
Consider a single particle in one dimension whose evolution is given by the NLSE, i ∂ψ ∂t = Hψ , (9.33) where (9.34) and p and v ext are the momentum of the particle and the external potential, respectively, and the dissipative potential is given by As can be easily verified, this NLSE has the very interesting property that it satisfies the Langevin equation at zero temperature for the expectation values of observables: Interestingly, v bath can be written as a functional of the particle density, This identification is very appealing, since the frictional force is proportional to the space integral of the velocity field of the particle, ĵ(z ′ ) t / n(z ′ ) t .Furthermore, the friction coefficient λ can be derived from a microscopic model of harmonic bath modes [Zwanzig 2001, Nitzan 2006, Tuckerman 2010].
Note that v bath at a given time only depends on the momentum of the particle at the same instant, implying that this NLSE is Markovian.This situation can be obtained in the limit where the dynamics of the bath can be described as white noise [Peskin 1998].Although the discussion above has been given for a single-particle, we heuristically propose Eq. ( 9.38) as a Markovian bath functional for TDDFT.In practice, we can re-express v bath in terms of the orbitals of a time-dependent single Slater determinant KS wavefunction, The extension of the functional to more dimensions follows analogously, although the limits of integration must be studied with care.Equation (9.39) is easy to implement in a real-time propagation, and has been implemented for a model Helium system interacting with a heat bath [Yuen-Zhou 2010].Non-Markovian extensions, as well as functionals where several timescales of relaxation and dephasing exist, are currently under development.
With slightly different motivations, Neuhauser and Lopata [Neuhauser 2008] have recently reported important result which could also be considered a Markovian bath functional in our formalism.Their functional is inspired by an optimal control approach, where they demand that the energy in the KS system decays monotonically.They show that achieves such goal.This functional couples the average current-density to the current operator with a spatially dependent proportionality constant a(z ′ ).
Their studies of a jellium cluster also show numerical robustness to provide dissipation in a real-time KS calculation.The merits and disadvantages of (9.38) and (9.40) shall be investigated in detail in the future.

OQS-TDDFT in the linear response regime using the open Kohn-Sham scheme
In addition to real-time dynamics, one can also consider OQS-TDDFT in the linear response regime, which gives access to environmentally broadened spectra [Tempel 2011a].The starting point is the density-density response function of an interacting OQS evolving according to Eq. (9.11): where ρn S (r, 0) = Tr R [n(r), ρeq ] (9.42) is the commutator of the particle density operator with the full equilibrium density matrix of the combined system and reservoir.From Eq. (9.41), one sees that the primary effect of the bath is to introduce a frequency-dependent self-energy K(ω), which shifts the poles of the response function into the complex plane.In the absence of coupling to the bath, the poles of Eq. (9.41) would lie at the excitation frequencies of the isolated system, which are the eigenvalues of LS .The real part of K(ω) can be interpreted as an excited-state lifetime while the imaginary part is a Lamb shift of the energy.
In order to access the poles of the many-body response function in Eq. ( 9.41), one introduces an auxiliary open, but non-interacting KS system with a density-density response function given by (9.43) Here, LKS is the Liouvillian for the ground or equilibrium-state Kohn-Sham-Mermin Hamiltonian [Kohn 1965] and ρKS S (0) is the corresponding KS density matrix.KKS (ω) is a KS self-energy which describes coupling of noninteracting electrons to the environment.It is chosen to be a one-body superoperator and easily constructed in terms of KS orbitals and eigenvalues of Lks .As mentioned in Sect.9.3, the existence and uniqueness of such a KS system is guaranteed by setting K′ = KKS , Ξ ′ = Ξ KS and v ′ ee (r, r ′ ) = 0 in the proof.For this scheme, we partition the potential as where v open xc not only accounts for electron-electron interaction within the system, but must also correct for the difference between KKS and K in the system-bath interaction.This scheme is better suited to response theory than the closed Kohn-Sham scheme discussed in Sect.9.4, since relaxation and dephasing is already accounted for in the KS system through KKS .The unknown (OQS-TDDFT) exchange-correlation functional only needs to correct the relaxation and dephasing in the KS system to that of the interacting system, rather than needing to explicitly account for the entire effect of the environment.However, the closed KS scheme is better suited for real-time dynamics, since one only needs to propagate a set of equations for the KS orbitals as in usual TDDFT.
Reminiscent of usual TDDFT, Eq. ( 9.41) and Eq. ( 9. is the OQS-TDDFT exchange correlation kernel.It is a functional of the equilibrium density as well as the memory kernel in both the interacting and Kohn-Sham systems.For Markovian environments, it is straightforward to reformulate Eq. ( 9.45) as a Casida-type equation [Casida 1996], (9.47)where the frequency-dependent operator Ω(ω) can be expressed in a basis of Kohn-Sham-Mermin orbitals as Here, Γ KS kl and ∆ KS kl arise from matrix elements of the real and imaginary parts of KKS respectively.
are matrix elements of the OQS Hartree-exchange-correlation kernel.Equation (9.48) is a non-hermitian and explicitly frequency-dependent operator yielding complex eigenvalues.The real part of the eigenvalues are interpreted as excitation energies while the imaginary parts give the linewidths.Since the KS system is open, the bare KS spectrum is already broadened at zerothorder.f open xc has the task of not only shifting the location of the bare KS absorption peaks to that of the interacting system, but it must also correct the linewidths.
As a simple example, we solved the OQS Casida equations in Eq. ( 9.47) to obtain the absorption spectrum of a C 2+ cation interacting with the modes of the electromagnetic field in vacuum, giving rise to radiative natural linewidths.The electromagnetic field acts as a photon bath, while the C 2+ cation can be treated as an OQS in our formalism [Cohen-Tannoudji 2004].As a crude first approximation, we used an adiabatic functional (ATDDFT) for f open xc in Eq. (9.48), and solved Eq. (9.47) for the 3 lowest dipole allowed transitions (2s → 2p, 3p, 4p).From Fig. (9.2), we see that the adiabatic functional places the location of the absorption peaks in essentially the correct place as in usual TDDFT, but leaves the linewidths unchanged relative to their bare KS value.To correct the linewidths, one needs a frequencydependent bath functional yielding additional broadening.Such a functional with the correct frequency-dependence to first-order in Görling-Levy perturbation theory [Görling 1993, Görling 1998a] was discussed in [Tempel 2011a] for the 2s → 2p transition.The frequency-dependent kernel matrix element in Eq. (9.49) was found to be where Γ KS 2p,2s is the bare KS linewidth and Γ 1 2p,2s is a correction derived from first-order Görling-Levy perturbation theory.In Fig. 9.3, we see that including Eq. (9.50) when solving the OQS Casida equations yields a large correction to The numerically exact spectrum obtained using experimental data.For visualization, all linewidths have been scaled by a factor of c 3 since the radiative lifetime in vacuum is extremely small.
the linewidth, although the oscillator strength is unchanged.To correct the oscillator strength as well, one needs higher-order corrections.A similar formalism can be used to capture line broadening due to vibrational relaxation in molecules and electron-phonon scattering and is currently being explored.

Conclusions and outlook
We have discussed the formal foundations of OQS-TDDFT in the density matrix representation starting from a many-electron quantum master equation and established a van Leeuwen construction which allows for a variety of different Kohn-Sham schemes.
The first scheme we discussed uses a non-interacting and closed (unitarily evolving) Kohn-Sham system to reproduce the dynamics of an interacting OQS.With suitable functionals, this scheme is remarkably useful for dissipative real-time dynamics, since it can be easily implemented in existing real-time codes.We presented the simple yet practical Markovian bath functional, which has shown promising results for a model Helium system [Yuen-Zhou 2010].Future research will focus on understanding exact conditions with the goal of developing more sophisticated functionals.In [Tempel 2011a], a systematic study of the exact OQS-TDDFT functional was carried out for a one-electron model OQS.The exact functional was shown to have memory dependence [Maitra 2001] and share some fea- tures with existing dissipation functionals in time-dependent current DFT (TDCDFT) [Vignale 1996, Vignale 1997, Ullrich 2002b].However, in OQS-TDDFT dissipation arises from coupling to a dense bosonic bath, which differs from TDCDFT where dissipation arises as an intrinsic feature of the interacting electron liquid.
The second scheme we discussed uses an open KS system to calculate broadened absorption spectra in linear response TDDFT.By using an open KS system, the bare KS spectrum is already broadened, while the OQS-TDDFT exchange-correlation kernel generates additional line broadening and shifts.The development of more sophisticated frequency-dependent functionals to capture additional broadening and asymmetric lineshapes due to non-Markovian effects is currently being explored as well.
In the next chapter, an alternative formulation of OQS-TDDFT based on stochastic wavefunctions rather than density matrices will be presented.
Fig. 9.3.The curves shown are: (a) The correction to the bare Kohn-Sham linewidth using the frequency-dependent bath functional to first-order in GL perturbation theory [Eq.(9.50)].(green).(b) The spectrum obtained by solving Eq. (9.47) with an adiabatic exchange-correlation kernel (blue-dashed).(c) The numerically exact spectrum obtained using experimental data (red-dashed).All linewidths are scaled by a factor of c 3