Time-Dependent Density Functional Theory of Open Quantum Systems in the Linear-Response Regime

Time-Dependent Density Functional Theory of Open Quantum Systems in the Linear-Response Regime

David G. Tempel Department of Physics, Harvard University, 17 Oxford Street, 02138, Cambridge, MA    Mark A. Watson    Roberto Olivares-Amaya    Alán Aspuru-Guzik Department of Chemistry and Chemical Biology, Harvard University, 12 Oxford Street, 02138, Cambridge, MA aspuru@chemistry.harvard.edu
Abstract

Time-Dependent Density Functional Theory (TDDFT) has recently been extended to describe many-body open quantum systems (OQS) evolving under non-unitary dynamics according to a quantum master equation. In the master equation approach, electronic excitation spectra are broadened and shifted due to relaxation and dephasing of the electronic degrees of freedom by the surrounding environment. In this paper, we develop a formulation of TDDFT linear-response theory (LR-TDDFT) for many-body electronic systems evolving under a master equation, yielding broadened excitation spectra. This is done by mapping an interacting open quantum system onto a non-interacting open Kohn-Sham system yielding the correct non-equilibrium density evolution. A pseudo-eigenvalue equation analogous to the Casida equations of usual LR-TDDFT is derived for the Redfield master equation, yielding complex energies and Lamb shifts. As a simple demonstration, we calculate the spectrum of a C atom in an optical resonator interacting with a bath of photons. The performance of an adiabatic exchange-correlation kernel is analyzed and a first-order frequency-dependent correction to the bare Kohn-Sham linewidth based on Görling-Levy perturbation theory is calculated.

pacs:

I I. Introduction

Due to its attractive balance between accuracy and efficiency, time-dependent density functional theory (TDDFT) has seen a tremendous growth of applications in recent years. These range from optical properties of molecules, clusters and solids, to optimal control theory and real-time dynamics of species in intense laser fields runge gross (); annu (); review of burke (); review (); rappoport (); li 1 (). TDDFT has been particularly successful at calculating optical response properties of electronic systems in the linear response regime linear response (). In most quantum chemical codes, excitation energies and oscillator strengths are extracted by solving a pseudo-eigenvalue equation, originally formulated by Casida  Casida Review (). The Casida equations are derived by considering the linear density response of an interacting system and corresponding non-interacting Kohn-Sham system, both undergoing unitary evolution. If the Casida equations are solved using the ubiquitous adiabatic approximation (ATDDFT) within a discrete basis set, the resulting eigenvalues are real. This gives rise to a discrete absorption spectrum of delta function peaks.

In experimentally observed spectra, line broadening arises from a variety of different mechanisms, several of which have been explored already within LR-TDDFT. In extended systems, relaxation and dephasing due to electron-electron scattering is well captured using non-adiabatic and current-density dependent functionals Vignale-Kohn (); VUC (); Vignale-Dgosta (); Ullrich 2-electron (). In finite systems, decay of resonant states due to coupling with the continuum gives rise to finite linewidths. The ability of DFT and TDDFT to capture lineshape parameters and widths of resonances has been discussed in neepa resonance (); exx resonance (); adam resonance ().

Another important broadening mechanism arises from relaxation and dephasing of electronic degrees of freedom by a classical or bosonic bath such as photons, phonons or impurities. For extended systems, this situation was considered in weakly disordered (); ullrich springer (); memory function (), where linewidths of intersubband plasmon excitations were well captured by combining the Vignale-Kohn functional Vignale-Kohn () to account for electron-electron scattering with the memory function formalism for electron-phonon and electron-impurity scattering. For atomic and molecular systems, the theory of open quantum systems within the master equation approach is often used nakajima (); zwanzig (); breuer (); Kossakowski72a (); lindblad (); gorini (). Several important examples include vibrational relaxation of molecules in liquids and solid impurities  Mukamel Non-Markov (); Mukamel Frank-Condon (); Nitzan () , cavity quantum electrodynamics (QED)  resonator (); spontaneous (); spontaneous1 (); spontaneous2 (), photo-absorption of chromophores in a protein bath Redfield Absorption (); Redfield Absorption 2 (); protien dynamics (); mukamel 1 (), single-molecule transport burke (); gebauer chemphyschem (); gebauer prl (); chen 1 (); chen 2 () and exciton transport exciton transfer (); nano letter (); ENAQT (); mukamel 2 (). In all of these examples, even with simple system-bath models, the exact solution of the master equation for the reduced dynamics of the many-body electronic system is computationally intractable. Therefore, open quantum systems TDDFT (OQS-TDDFT) offers an attractive approach to the many-body open-systems problem.

Several important steps toward the formulation of OQS-TDDFT have recently been made, with the focus on real-time dynamics. In burke (), a Runge-Gross theorem was established for Markovian master equations of the Lindblad form. A scheme in which the many-body master equation is mapped onto a non-interacting Kohn-Sham master equation was proposed for application to single-molecule transport. In joel's paper (); prl (), the Runge-Gross theorem was extended to arbitrary non-Markovian master equations and a Van Leeuwen construction was established, thereby proving the existence of an OQS-TDDFT Kohn-Sham scheme van leeuwen (). In prl (), it was shown that the original interacting open system dynamics can be mapped onto either a non-interacting open Kohn-Sham system, or a non-interacting closed (unitarily evolving) Kohn-Sham system. A different formulation of OQS-TDDFT based on the stochastic Schrodinger equation has also been developed d prl (); da prb (); d prb ().

The goal of the present manuscript is to formulate the linear response version of (OQS-TDDFT) within the master equation approach. This provides a framework in which environmentally broadened spectra of many-body electronic systems can be accessed in an ab initio way using TDDFT, especially when combined with microscopically derived master equations. We use the scheme discussed in prl (); burke (), in which the interacting OQS can be mapped onto a non-interacting open Kohn-Sham system, yielding the same density response. This scheme is better suited to response theory than the closed Kohn-Sham scheme also discussed in prl (), since relaxation and dephasing is already accounted for in the Kohn-Sham system. The unknown (OQS-TDDFT) exchange-correlation functional only needs to correct the relaxation and dephasing in the Kohn-Sham system to that of the interacting system, rather than needing to explicitly account for the entire effect of the environment. However, the closed Kohn-Sham scheme is better suited for real-time dynamics, since one only needs to propagate a set of equations for the Kohn-Sham orbitals as in usual TDDFT. This is in contrast to the open Kohn-Sham scheme, in which equations are propagated for the elements of the density matrix, with N being the dimensionality of the Hilbert space.

The paper is organized as follows. In Section II, we formulate the most general OQS-TDDFT linear response equations for arbitrary non-Markovian master equations with initial correlations. In section III, we make the treatment more specific by focusing on the Redfield master equation. We also derive Casida-type equations whose solution yields the environmentally broadened absorption spectrum. The solutions to these equations are complex, with the real part of the frequency yielding the location of absorption peaks and the imaginary part yielding the linewidths. In Section IV, we apply the formalism developed in Section III to a C atom in an optical resonator setup evolving under the Redfield master equation. Section V analyzes the performance of using an adiabatic functional (OQS-ATDDFT) in solving the OQS-TDDFT Casida equations derived in section III. To a large degree, OQS-ATDDFT is seen to provide a reliable correction to the location of absorption peaks while leaving the linewidths unchanged. A frequency-dependent functional yielding a correction to the OQS-ATDDFT linewidth based on Görling-Levy (GL) perturbation theory is then calculated and analyzed. In section VI a discussion and outlook is provided.

Atomic units in which are used throughout. This also implies that the speed of light in vacuum is given by . For generality, we have formulated most of the theory by considering linear response from an equilibrium state at finite temperature. For atoms and molecules, it will generally be sufficient to take the zero temperature limit and consider linear response from the ground-state.

Ii II. General formulation of OQS-TDDFT linear response theory

A. Linear response of interacting open quantum systems

Our formulation of the interacting OQS density-density response function parallels that used in spin response () for calculating spin susceptibilities (see also Linear Response Open (); reduced response ()). The starting point is the unitary evolution for the full density matrix of the system and the reservoir (we use the terms "reservoir" and "bath" interchangeably throughout) ,

(1)

where is the Liouvillian superoperator for the full system and reservoir dynamics. The full Hamiltonian is given by

(2)

Here,

(3)

is the Hamiltonian of the electronic system of interest in an external potential . 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. The system-bath coupling, , is at this point arbitrary, but we will discuss specific forms later. acts in the combined Hilbert space of the system and reservoir and so it couples the two subsystems. is the Hamiltonian of the reservoir, assumed to have a dense spectrum of eigenstates. The density of states of determines the structure of reservoir correlation functions, whose time-scale in turn determines the reduced system dynamics.

Defining the reduced density operator for the electronic system alone by tracing over the reservoir degrees of freedom,

(4)

one arrives at the formally exact quantum master equation,

(5)

Here, is the memory kernel and arises from initial correlations between the system and its environment. It is referred to as the inhomogeneous term. The above equation is still formally exact, as gives the exact expectation value of any observable depending only on the electronic degrees of freedom. In practice, however, approximations to and are required. Of particular importance in TDDFT is the time-dependent electronic density,

(6)

where . We now assume that for , the external potential is time-independent while for a weak perturbing field is applied. i.e.

(7)
(8)

For , the entire system and environment is in thermal equilibrium described by the canonical density operator

(9)

where is the inverse temperature. The reduced equilibrium density operator of the electronic system is then given by

(10)

In Eq. 9 and Eq. 10, is the full Hamiltonian for and is the static Hamiltonian of the electrons in the absence of the external perturbation.

For , the perturbing field is switched on and the system density operator subsequently evolves under the master equation given in Eq. 5. The electronic density evolution to first-order in the perturbing field is then given by

(11)
(12)

Here, is the equilibrium electron density and

(13)

is the linear density response. is the density-density response function. Its Fourier transform to the frequency domain is given by

(14)

where

(15)

is the operator generating the electronic charge density in the Heisenberg picture with respect to the full Hamiltonian for . Rearranging terms under the trace operation in Eq. 14, the density-density response function can be written as

(16)

Here, is an operator acting in the electronic system Hilbert space, which obeys the same equation of motion as the reduced system density operator (Eq. 5)

(17)

subject to the initial condition

(18)

Carrying out the Fourier transform in Eq. 16, one arrives at the formally exact expression for the open-systems density-density response function in Liouville space,

(19)

is the Liouvillian for the system Hamiltonian for , defined by it’s action on an arbitrary operator by . It is readily verified that Eq. 19 reduces to the usual expression for the density-density response function of an isolated system when , and .

The absorption spectrum can be extracted by taking the imaginary part of in Eq. 19. For an isolated system with a discrete spectrum, this is given by a sum over weighted delta function peaks. For an open system as in Eq. 19, and in principle give rise to the exact complicated broadened and shifted spectrum, due to relaxation and dephasing of the electronic degrees of freedom by the environment. In practice, however, even with simple approximations to and , the exact form of is not exactly known, since it refers to a many-body response function. In the next two subsections, we consider an open Kohn-Sham system, formally yielding the exact density response of the original interacting open system. In an open-systems TDDFT framework, the exact spectrum of Eq. 19 is obtained by correcting the open Kohn-Sham spectrum via an exchange-correlation kernel. The kernel must take into account not only the electron-electron interaction contained in , but must also correct the interaction of the system with the bath, described by and .

B. The open Kohn-Sham system

It was proven in prl (), that for a master equation of the form given in Eq. 5, there exists a unique, non-interacting and open Kohn-Sham system, whose system density operator evolves under the master equation

(20)

such that the time-dependent density is obtained from

(21)

for all times. , where the Kohn-Sham Hamiltonian is given by

(22)

Here, is a local, multiplicative, one-body potential which drives the open Kohn-Sham system in such a way that the true density of the original interacting open system is reproduced for all times. In analogy to usual TDDFT, the Kohn-Sham potential is partitioned as

(23)

where is the Hartree potential and the unknown functional accounts for electron-electron interaction within the system as well as interaction between the system and bath. In general,

(24)

Formally, the open-systems exchange-correlation potential is a functional not only of the density, but also of the memory kernel, inhomogeneous term and initial state of both the interacting and Kohn-Sham systems. It has been shown in usual TDDFT of closed systems, that initial state dependence can be absorbed as dependence on the history of the density and vice versa neepa history (). Interestingly, in the theory of open quantum systems, it is possible to absorb the inhomogeneous term , into the memory kernel  tannor (). This raises the possibility that may be a functional only of n, and , but a more rigorous study of this will be done in future work. For notational convenience, we suppress the explicit functional dependence of on these quantities, although it is implied unless otherwise stated.

In general, and can be chosen to simplify as much as possible, although with some restrictions prl () and consistency conditions between and tannor ().

If the system is started in an equilibrium state, as is typically the case in linear response theory, the initial state dependence in Eq. 24 is automatically removed. The equilibrium density, , is obtained by solving the Kohn-Sham-Mermin equations kohn-sham ()

(25)

The Kohn-Sham-Mermin potential is partitioned as

(26)

where is the exchange-correlation contribution to the free energy. After solving Eq. 25, the equilibrium Kohn-Sham-Mermin density operator is obtained by populating the orbitals according to

(27)

where are Fermi-Dirac occupation numbers

(28)

Denoting , the one-particle Kohn-Sham-Mermin density matrix is

(29)

The equilibrium density is then obtained by taking the diagonal elements in real space,

(30)

For the open Kohn-Sham scheme to be useful practically, , and should be constructed so that the following conditions are satisfied:

1) and should not induce correlations between non-interacting electrons as the Kohn-Sham system evolves. This ensures that the N-body Kohn-Sham density matrix in Eq. 20 can be traced over N-1 electron coordinates to arrive at a closed equation of motion for the Kohn-Sham reduced 1-particle density matrix. Physically, this is expected since most reasonable bath models couple to the electronic system through one-body operators Cohen (); breuer ().

2) The equation of motion for the non-equilibrium open Kohn-Sham reduced 1-particle density matrix should have the Kohn-Sham-Mermin density matrix as its stationary-point solution. This ensures that the system thermalizes to the correct equilibrium density. This also means that at equilibrium, should reduce to . Although this is automatically satisfied by using as an adiabatic approximation, it might not be satisfied by more sophisticated approximations with memory dependence neepa history (); neepa memory (); neepa memory perturbations (); baer memory ().

C. Linear response of the open Kohn-Sham system

Returning to linear response, for the Kohn-Sham system is in thermal equilibrium with its environment at inverse temperature , described by Eq. 25. At , the perturbation is switched on and the Kohn-Sham system subsequently evolves according to Eq. 20. The Kohn-Sham potential is given by

(31)
(32)

where . Due to Eq. (21), the exact linear density response of Eq. (12) is obtained through

(33)

Here, is the density-density response function of the open Kohn-Sham system. It’s Fourier transform to the frequency domain is given by

(34)

where is the Liouvillian for the equilibrium Kohn-Sham-Mermin Hamiltonian. Since the system is in the equilibrium state at , must yield the equilibrium density, implying that it reduces to the Kohn-Sham-Mermin density matrix when traced over N-1 electron coordinates. We now define the open-systems exchange-correlation kernel in analogy to usual TDDFT for closed-systems by,

(35)

which is a functional of the equilibrium density. As in TDDFT for closed systems, the interacting and Kohn-Sham response functions are related through a Dyson-like equation,

(36)

, being much simpler then the original interacting , can readily be constructed from the orbitals and eigenvalues in Eq. 25 and approximations to and in terms of these quantities. This will be done explicitly for the Redfield master equation in the next section. Since correlation between the open Kohn-Sham system and reservoir is already partially captured through and , the bare Kohn-Sham absorption spectrum extracted from is already broadened and shifted. The functional has the task of correcting the spectrum extracted from to that of the interacting , incorporating both the usual electron-electron correlation in closed-systems TDDFT as well as additional system-bath correlation. In general, the memory kernel may give rise to a very complicated non-analytic structure of in the lower half of the complex plane. However, for Markovian master equations, it will be seen that the pole structure of in the discrete part of the spectrum consists of simple poles in the lower complex plane, shifted by a finite amount off of the real axis. In such cases, it might be reasonable to approximate as

(37)

Here, the first term is just the adiabatic contribution to the exchange-correlation kernel and the second term is an in general frequency-dependent and imaginary correction. This is attractive, since we can take advantage of the usual good performance of adiabatic TDDFT in describing the location of absorption peaks, and attempt to build functionals that go beyond the adiabatic approximation to account for line broadening and lamb shifts. This strategy will be discussed further in section V.

Iii III. LR-TDDFT for the redfield master equation.

In section II, we formulated LR-TDDFT for a very general class of master equations. In this section, we make the discussion more specific by invoking the Markov approximation and second Born approximation in the system-bath interaction, to arrive at the microscopically-derived Redfield master equation Redfield Absorption (); Cohen (); May/Kuhn (); van Kampen (); redfield original 2 (); redfield original (). Since the Redfield equations are rigorously obtained without phenomenological parameters, they are amenable to an ab initio theory such as TDDFT. Although we focus on Redfield theory here, the generalization of our formalism to other Markovian master equations can be made with small modifications. Finally, we discuss how it is possible to extract the absorption spectrum of a many-body system evolving under the Redfield equations directly within OQS-TDDFT. This is done by formulating Casida-type equations yielding complex eigenvalues due to coupling with the bath.

A. The Markov approximation and the Redfield master equation.

The Markov approximation describes a situation in which the bath correlation functions decay on an infinitely fast time-scale relative to the thermalization time of the system May/Kuhn (); Cohen (). As a result, the bath has no memory and the memory kernel is time-local

(38)

Additionally, this implies that the initial density operator is a tensor product of a density operator in the system space with the equilibrium density operator of the bath

(39)

As a result of Eq. 39, the system and environment have no initial correlations, and

(40)

The master equation (Eq. 5) then takes the simple form,

(41)

If the system Hamiltonian is time-independent, Eq. 41 is written in a basis of eigenstates of as:

(42)

Here, are many-body transition frequencies of and are matrix elements of in this basis. So far our discussion applies to any Markovian master equation. To obtain the Redfield equations, we further assume that the system-bath coupling has a bilinear form

(43)

where is an operator in the reservoir Hilbert space which couples to a local one-body operator in the system Hilbert space. This form of the system-bath coupling is very general and can apply to electron-phonon coupling in molecules and solid impurities, but also momentum dependent couplings which are relevant for instance in laser cooling, brownian motion in liquids or dissipative strong field dynamics tannor (); Kohen (); Cohen (). The Redfield tensor is then derived by performing second-order perturbation theory in , and is given explicitly by

(44)

For a detailed derivation of the Redfield equations see May/Kuhn (). In Eq. 44,

(45)

are matrix elements of between system many-body wavefunctions and

(46)

are system many-body excitation energies. are bath correlation functions given by

(47)

where

(48)

B. Linear response of a many-body system evolving under the Redfield master equation.

Since we consider linear response from the equilibrium state, the initial density matrix for the system is given by

(49)

and the density-density response function in Eq. 19 reduces to

(50)

Inserting a complete set of eigenstates of in Eq. 50, a sum over states expression for the density-density response function is given by,

(51)

Here, are equilibrium occupation probabilities of the various many-body states.

By hermiticity of the density matrix, it can be readily verified that . We can separate the real and imaginary parts of as

(52)

From the pole structure of Eq. 51, we see that corresponds to an imaginary part of the energy of the transition , giving rise to a finite lifetime, while is a Lamb shift of the real part of the energy. The effect of the Redfield tensor is to shift the poles of the density-density response function by a finite amount into the lower half of the complex plane.

C. The Markovian Kohn-Sham-Redfield equations

We now discuss the properties of the open Kohn-Sham system for the Redfield master equation. As discussed in section II B, there is some freedom in the construction of the Kohn-Sham dissipative superoperator and corresponding exchange-correlation potential. In this section, we choose a very natural form for the Kohn-Sham superoperator, whose construction is discussed below.

We consider an open Kohn-Sham system evolving under a Markovian master equation

(53)

which reproduces the exact density evolution of the interacting Markovian master equation in Eq. 41. To satisfy condition 1 in section II B, we must have

(54)

i.e. the N-body Kohn-Sham dissipative superoperator is a sum of one-body superoperators acting on each coordinate separately. This also follows very naturally from the assumed one-body nature of the system-bath interaction in Eq. 43. Since is also a sum of one body terms, we can trace both sides of Eq. 53 over N-1 electron coordinates and arrive at a closed equation of motion for the Kohn-Sham 1-particle reduced density matrix,

(55)

We can now write Eq. 55 in a basis of Kohn-Sham-Mermin orbitals as

(56)

We choose to have the form of the Redfield tensor, but written in terms of Kohn-Sham-Mermin orbitals and eigenvalues,

(57)

Here,

(58)

are bare Kohn-Sham-Mermin transition frequencies and

(59)

are matrix elements of the system-bath coupling operator between Kohn-Sham-Mermin orbitals.

Eq. 56 has a number of desirable properties. First, it’s stationary point solution is the Kohn-Sham-Mermin density matrix, Eq. 27, if reduces to the Kohn-Sham-Mermin Hamiltonian when evaluated on the equilibrium density. This ensures that condition 2 in section II B is satisfied. It also satisfies detailed balance as well as most other properties of the usual many-body Redfield equations, but in terms of Kohn-Sham-Mermin quantities. Also, the tensor has a simple form and can be constructed explicitly in terms of orbitals and eigenvalues obtained in an equilibrium-state Kohn-Sham calculation. The potential contained in will in general be a functional of and as well as the time evolving density.

D. Linear response of the Kohn-Sham-Redfield system and the open-systems Casida equations

We now consider the open-systems LR-TDDFT formalism developed in section II, but applied to the Redfield master equation. The density-density response function of the Kohn-Sham-Redfield system is given by

(60)

Using Eqs. 53-56 and inserting a complete set of Kohn-Sham-Mermin states, one obtaines the sum-over-states expression,

(61)

Eq. 51 and Eq. 61 are related through the Dyson-like relation given in Eq. 36. To extract the poles of the interacting density-density response function in Eq. 51 from that of the Kohn-Sham system in Eq. 61, a pseudo-eigenvalue equation must be solved for the squares of the complex transition frequencies,

(62)

The operator can be written as a matrix in a basis of Kohn-Sham molecular orbitals (assuming a closed shell system) as

(63)

The explicit derivation of Eq. 62 and Eq. 63 is given in appendix A. In Eq. 63,

(64)

and the bare Kohn-Sham linewidths and Lamb shifts are given by the relation

(65)

as in Eq 52. In principle, with the exact functional , the exact poles of Eq. 51 are recovered by solving Eq. 62. The operator is non-Hermitian, giving rise to complex eigenvalues corresponding to broadened excitation spectra. is also frequency-dependent and imaginary, both explicitly through the third term in Eq. 63 and implicitly through in the coupling matrix . The explicit frequency-dependence arises because the bare Kohn-Sham transitions are already broadened, even in the absence of Hartree-exchange-correlation effects. This is most easily seen by setting in Eq. 63. is then diagonal and Eq. 62 reduces to a set of uncoupled equations given by

(66)

These are solved with the quadratic formula to yield

(67)

which are precisely the poles of Eq. 61.

Iv IV. Application - Spectrum of a C atom from the Redfield master equation

As a simple demonstration, in this section we calculate the absorption spectrum of an atom in an optical resonator using the Redfield master equation Cohen (); breuer (). The modes of the radiation field act as a bosonic bath, leading to decay of the atomic excitations due to spontaneous and stimulated emission resonator (); spontaneous1 (); spontaneous2 (); spontaneous (). We consider a hypothetical experimental setup similar to that used in resonator (). The effect of the optical resonator is to modify the density of radiation field modes relative to the vacuum, leading to an enhancement of the decay rate of atomic excitations at certain frequencies resonator (); Cohen (); quantum optics review (). Although we do not have access to experimental data for C in an optical resonator cavity, we use accurate experimental data taken in vacuum, together with a chosen cavity geometry to construct a numerically exact spectrum NIST ().

Spontaneous emission from the Redfield master equation

For a single atom with zero center of mass velocity contained in an optical resonator cavity, the system-bath interaction is

(68)

Here, is the dipole operator for the atom, and are respectively the polarization vector and frequency of the ith mode of the radiation field and V is the quantization volume. and respectively destroy and create a photon in the ith mode of the cavity. The photon reservoir Hamiltonian is

(69)

With the system-bath interaction and bath hamiltonian specified, the Redfield tensor can be explicitly constructed Cohen ().

The entire atom-field system is taken to be in thermal equilibrium at inverse temperature , such that . With this condition, the atom can be assumed to be in it’s groundstate and the effect of stimulated emission is neglected. We then need to only construct the matrix elements appearing in Eq. 51. For the real part of one finds

(70)

The imaginary part of is given by

(71)

where denotes the principle value integral Cohen (). is the frequency of a photon whose wave-vector magnitude is . is the difference in atomic energy levels and are matrix elements of the dipole operator between atomic wavefunctions. is the density of field modes in the cavity. In free-space, the density of field modes takes the form . This gives rise to a free-space natural linewidth

(72)

If one considers an experimental setup such as that used in resonator (), the density of field modes is modified to

(73)

where the function is given by

(74)

The linewidth in the cavity is then modified to

(75)

where we have neglected cavity edge effects and angular dependence of . In Eq. 74, L is the length of the cavity and , where R is the reflection coefficient of the cavity walls. By changing the mirror reflectivity and cavity length, suppression or enhancement of the spontaneous emission rate is possible. For our calculations, we choose cavity parameters of and cm, leading to an overall enhancement. Using experimental data for the atomic energy levels and the transition dipole matrix elements of C taken from NIST (), together with the specified parameters for the cavity geometry, we can explicitly construct Eq.  75. We include in our calculation the 3 lowest dipole allowed transitions of C, which are . The numerical values of the linewidths (imaginary part of the frequency) calculated in Eq. 75 are given in the fourth column of Table 2. Due to lack of experimental data on all atomic transitions, the Lamb shifts in Eq. 71 cannot be explicitly evaluated. They are, however, estimated to be several orders of magnitude smaller and will be neglected in the following analysis. From the transition dipole matrix elements we can also construct the oscillator strengths. The "exact" response function constructed with these parameters is included in Figures 1-3.

OQS-ATDDFT calculation of the spectrum of C

As a first step, we solve Eq. 62 using only an adiabatic approximation to the exchange-correlation kernel. This corresponds to including the first term in Eq. 37, but entirely neglecting .

For our calculations, we obtain the Kohn-Sham parameters of C to be inputted in Eq. 63 using the real-space TDDFT package Octopus octopus 1 (); octopus 2 (); octopus 3 (). First, a ground-state DFT calculation is performed using the local density approximation (LDA) with the modified Perdew-Zunger (PZ) parameterization of the correlation energy PZ (). For all calculations, the core is replaced by a Troullier-Martins pseudopotential TM (). The Kohn-Sham eigenvalues and dipole matrix elements between Kohn-Sham orbitals are computed, and substituted into Eq. 75 to obtain the bare Kohn-Sham linewidths, . These correspond to the real part of the Kohn-Sham-Redfield tensor (imaginary part of the bare Kohn-Sham frequency) of Eq. 57 and are given in the second column of Table 2. From the dipole matrix elements between Kohn-Sham orbitals, the bare Kohn-Sham oscillator strengths can be constructed. The real and imaginary parts of the bare Kohn-Sham response function constructed with these quantities are plotted in figures 1-3.

Next, we perform a standard LR-ATDDFT calculation to obtain the matrix elements of the adiabatic Hartree-exchange-correlation kernel to be inputted in Eq. 63. The matrix elements of the kernel obtained are given in table 3. We also include the energies obtained from the standard LR-ATDDFT calculation in column 4 of Table 1. For consistency, we use the LDA with modified PZ functional for the exchange-correlation kernel as well.

Since the operator in Eq. 63 is explicitly frequency-dependent even when using an adiabatic kernel, Eq. 62 represents a non-linear eigenvalue problem. We solve it using the generalized eigenvalue algorithm presented in generalized 1 (); generalized 2 (). The real part of the solutions to Eq. 62 are given in column 3 of Table 1, while the imaginary part is given in column 3 of Table 2. The real and imaginary parts of the response function obtained are plotted in figures 1-3.

Figure 1: Absorption Spectrum of C including the 3 lowest dipole allowed transitions. The curves shown are: a) The bare Kohn-Sham spectrum (green-dashed). b) The spectrum obtained by solving Eq. 62 with an adiabatic exchange-correlation kernel (blue-solid). c) The numerically exact spectrum obtained using experimental data (red-dashed). Also shown is the stick spectrum obtained by solving the usual Casida equations for C using ALDA (black-dotted).
Figure 2: Same as Figure 1, but with a close-up view of the and transitions.
Figure 3: Real part of the density-density response function of C including the 3 lowest dipole allowed transitions. The effect of the photon bath is to broaden the dispersion over multiple frequencies. The curves shown are: a) The bare Kohn-Sham dispersion (green-dashed). b) The dispersion relation obtained by solving Eq.62 with an adiabatic exchange-correlation kernel (blue-solid). c) The dispersion relation obtained using experimental data (red-dashed).
Transition Bare Kohn-Sham OQS-ATDDFT ATDDFT Exact
2s 2p 0.311 0.443 0.443 0.467
2s 3p 1.116 1.107 1.107 1.180
2s 4p 1.368 1.361 1.363 1.470
Table 1: Real part of the 3 lowest transition frequencies for C in an optical resonator in a.u.
Transition Bare Kohn-Sham OQS-ATDDFT Exact OQS-ATDDFT + GL
2s 2p 1.329 1.331 2.932 1.805
2s 3p 5.162 5.159 5.740
2s 4p 2.443 2.444 1.089
Table 2: Imaginary part of the 3 lowest transition frequencies for C in an optical resonator in a.u. The last column includes the GL perturbation correction to the 2s 2p transition.

V V. Analysis

Effect of using an adiabatic approximation to

From Tables 1 and 2, it is clear that using an adiabatic kernel in Eq. 62 gives essentially the same corrections to the real part of the energy as usual LR-ATDDFT, while giving almost no correction to the imaginary part. This means that using an adiabatic kernel in OQS-TDDFT is expected to give the same reliable correction to the location of absorption peaks as usual LR-ATDDFT, while correcting the bare Kohn-Sham linewidths requires an additional frequency-dependent functional. This justifies a posteriori our separation of the exchange-correlation kernel in Eq. 37 into an adiabatic part and a frequency-dependent part due exclusively to bath effects.

To understand this situation better, we consider a "small matrix approximation" (SMA), in which a single occupied to unoccupied Kohn-Sham transition, , is completely isolated from all other transitions Appel SMA (); Ullrich 2-level (); Neepa SMA (); Appel DPA (). This is valid when the transition of interest is weakly coupled to all other excitations. In this case, Eq. 62 reduces to