Fourth order expressions for the electronic absorption lineshape of molecular excitons

Jang, Seogjoo J.

doi:10.1063/1.5100986

The line shape of electronic absorption spectroscopy reflects the information on quantum dynamical processes accompanying the electronic excitation, and its accurate description is an important component for validating theoretical models and assumptions. The present work provides detailed expressions for the absorption line shape of molecular excitons that are valid up to the fourth order of exciton-bath interactions within the quantum master equation approach. These expressions can serve as the basis for developing general and systematic methods to model the line shape for a broad class of molecular exciton systems and environments. For the bath model of linearly coupled harmonic oscillators, more detailed expressions employing the spectral densities of the bath are presented. These expressions are then tested for a linear aggregate of identical chromophores each coupled to the harmonic oscillator bath. Calculation results for a super-Ohmic spectral density with exponential cutoff demonstrate the feasibility of calculations and also offer insights into the utility and difficulty of going beyond the second order approximation.

I. INTRODUCTION

Electronic absorption spectroscopy, one of the simplest experimental probes of quantum dynamical processes accompanying an electronic excitation conveys rich dynamical information through its line shape in the frequency domain.^1,2 However, compared to its routine experimental practice, accurate theoretical calculation of the line shape of an electronic absorption spectroscopy still remains difficult, in particular, for large molecular systems in complex environments. This is because the theoretical challenge entailed in such a calculation is the one commonly encountered for general open system quantum dynamics calculation in condensed media.^2–5

The analysis of the absorption line shape has played a central role in the spectroscopy of molecular excitons.^2,5–14 For well ordered molecular aggregates with relatively weak electron-phonon couplings, certain features of the line shape have played a central role in identifying unique features of the collective quantum nature of molecular excitons.^6,9,14 On the other hand, in disordered/fluctuating environments and in the presence of significant electron-phonon coupling, the mixing of inhomogeneous and homogeneous broadening mechanisms tends to blur such features, rendering unambiguous identification of excitonic peaks much more difficult.⁵ This has motivated the development of nonlinear (NL)^2,15,16 and single molecule (SM) spectroscopic techniques,^17–20 which can offer more selective information.

Despite great advances, analyses of NL or SM spectroscopic techniques have mostly relied on either simple models or additional set of assumptions that can introduce extra sources of ambiguity. On the contrary, the linear absorption spectroscopy, being the most universal and unambiguous experimental probe, requires the least amount of assumptions for its interpretation and thus continues serving as an essential tool for validating theoretical models and assumptions. In addition, line shape functions can also serve as important means for modeling exciton dynamics within the multichromophoric Förster resonance energy transfer theory²¹ or its non-Markovian generalization.²²

During the past decade, there have been significant advances in theories and computational methods for quantum dynamical calculation of spectral line shapes. The quantum master equation (QME) approach, from which the line shape can be obtained as the average of a dipole-dipole correlation function, has played a major role in this development.^{3–5,23–29} QMEs at the 2nd order level are now fairly well established^4,5,30–34 and can be applied to a wide range of systems. Significant advances have also been made in higher order approximations^35–38 and nonperturbative approaches,^26,38–45 although they have been limited to rather simple models and conditions.

For general models and environments, there is still great demand for further development of the QME approach beyond the 2nd order level. For accurate modeling of large systems in disordered environments, theories and computational methods need to be both accurate and efficient. In particular, an efficient method is important for disordered systems, where reliable sampling over the ensemble of disorder is necessary. An important avenue to explore in achieving this objective is systematic improvement of standard perturbative approximations with the QME approach,^35–37,46 at both theoretical and computational levels. Outcomes of these can also be utilized for further improvement through extrapolation schemes.^45,47–50

The main objective of the present work is to provide general expressions for the absorption line shape of molecular excitons that are valid up to the 4th order of the exciton-bath interactions. Particular attention is paid to deriving a general expression that can be directly employed for practical calculations. Model calculations are provided, which demonstrate the feasibility of applying the derived expressions. These model calculation results offer some promise but also indicate the difficulty of truncating expressions at a higher order without appropriate closure.

The paper is organized as follows. Section II provides the formalism in a manner more general than conventional ones and includes the effects of collective time dependent fluctuations of ground electronic state and average exciton energies. More detailed expressions are presented in Sec. III, and model calculations for a linear chain of identical chromophores, each coupled to an independent bath of harmonic oscillators, are provided in Sec. IV. Conclusions are provided in Sec. V.

II. GENERAL FORMALISM

A general formalism of the absorption line shape is presented. This is somewhat different from the more widely known dipole-dipole correlation function approach. A portion of this formulation has been presented⁴⁶ before in different notations, which is now incorporated into a more comprehensive formalism as described below.

Let us consider the combination of chromophores and host media, which are called material and represented by the following exciton-bath Hamiltonian:

Ĥ_{m} (t) = E_{g} (t) | g (⟩ ⟨) g | + \sum_{j} E_{e, 0} (t) | s_{j} (⟩ ⟨) s_{j} | + Ĥ_{e} + Ĥ_{e b} + Ĥ_{b},

(1)

where |g⟩ is the ground electronic state and E_g(t) is its electronic energy. This is assumed to be time dependent in general due to the influence of external agents or environmental effects that are not accounted for by the bath Hamiltonian Ĥ_b. This added time dependence can represent fluctuations due to local electric fields or effects of distant environments such as membranes for light harvesting complexes,⁵ which are not representable by simple bath Hamiltonians. In the second term, |s_j⟩ is the excitation state of a chromophore localized at site j, and E_e,0(t) is the average of the energies of |s_j⟩’s. The time dependence of E_e,0(t) has similar sources as that for E_g(t). Namely, they can result from external agents or environmental effects unaccounted for by the bath Hamiltonian. Spontaneous decays of |s_j⟩’s are assumed to take much longer time and are not included here.

The third term Ĥ_e in Eq. (1) is the exciton Hamiltonian representing relative energies and couplings of chromophores and is assumed to have the following form:

Ĥ_{e} = \sum_{j = 1}^{N_{c}} E_{j} | s_{j} (⟩ ⟨) s_{j} | + \sum_{j = 1}^{N_{c}} \sum_{k = j + 1}^{N_{c}} J_{j k} (| s_{j} (⟩ ⟨) s_{k} | + | s_{k} (⟩ ⟨) s_{j} |) .

(2)

In the above expression, E_j is the relative energy [with respect to E_e,0(t)] of the state |s_j⟩, and J_jk is the electronic coupling between states |s_j⟩ and |s_k⟩. It is assumed that |s_j⟩’s above are independent of nuclear and bath degrees of freedom (being based on the crude adiabatic approximation). It is also assumed that the variation of the relative energies of these states and their electronic couplings can be solely accounted for by the following exciton-bath Hamiltonian:

Ĥ_{e b} = \sum_{j = 1}^{N_{c}} \sum_{k = 1}^{N_{c}} {\hat{B}}_{j k} | s_{j} (⟩ ⟨) s_{k} |,

(3)

where ${\hat{B}}_{j k}$ is the bath component of the coupling to the electronic term |s_j⟩⟨s_k| and satisfies the condition ${\hat{B}}_{j k} = {\hat{B}}_{k j}^{†}$ ⁠.

For the material Hamiltonian given by Eq. (1), the time evolution operator can be expressed as

\begin{matrix} Û_{m} (t, t^{'}) & = & \exp_{(+)} \{- \frac{i}{ℏ} \int_{t^{'}}^{t} d τ Ĥ_{m} (τ)\} \\ = & | g (⟩ ⟨) g | U_{g} (t, t^{'}) e^{- i Ĥ_{b} (t - t^{'}) / ℏ} + U_{e, 0} (t, t^{'}) {\hat{1}}_{e} Û_{e b} (t - t^{'}) {\hat{1}}_{e}, \end{matrix}

(4)

where the subscript (+) represents chronological time ordering of the exponential operator and ${\hat{1}}_{e} = \sum_{j} | s_{j} (⟩ ⟨) s_{j} |$ is the unit operator defined in the excitonic subspace of electronic degrees of freedom. Time evolution functions and operator introduced above are defined as follows:

U_{g} (t, t^{'}) = \exp \{- \frac{i}{ℏ} \int_{t^{'}}^{t} d τ E_{g} (τ)\},

(5)

U_{e, 0} (t, t^{'}) = \exp \{- \frac{i}{ℏ} \int_{t^{'}}^{t} d τ E_{e, 0} (τ)\},

(6)

Û_{e b} (t) = \exp \{- \frac{i}{ℏ} (Ĥ_{e} + Ĥ_{e b} + Ĥ_{b}) t\} .

(7)

In addition, as will become clear below, it is convenient to introduce the following time evolution operator for the exciton Hamiltonian:

Û_{e} (t) = \exp \{- \frac{i}{ℏ} Ĥ_{e} t\} .

(8)

In the presence of the external radiation described at the semiclassical level, the total Hamiltonian can be expressed as

Ĥ_{T} = Ĥ_{m} (t) + Ĥ_{e r} (t),

(9)

where, for plane waves with electric polarization e and angular frequency ω_r, the exciton-radiation interaction Hamiltonian Ĥ_er(t) can be expressed as

Ĥ_{e r} (t) = \sum_{j = 1}^{N_{c}} \{A (μ_{j} \cdot e) e^{- i ω_{r} t} | s_{j} (⟩ ⟨) g |) (+ A^{*} (μ_{j} \cdot e^{*}) e^{i ω_{r} t} | g (⟩ ⟨) s_{j} |\} .

(10)

The above expression is based on the assumption that spatial variation of radiation at the length scale of the molecular system is negligible. Introducing

| D_{e} (⟩ = \sum_{j = 1}^{N_{c}} (μ_{j} \cdot e) | s_{j} (⟩,

(11)

which is the sum of transition dipole weighted excitation states, Eq. (10) can also be expressed as

Ĥ_{e r} (t) = \sum_{j = 1}^{N_{c}} \{A e^{- i ω_{r} t} | D_{e} (⟩ ⟨) g | + A^{*} e^{i ω_{r} t} | g (⟩ ⟨) D_{e} |\} .

(12)

For t ≤ 0, it is assumed that there is no radiation and that the material is in thermal equilibrium at temperature T for the ground electronic state. Then, the density operator representing the material at time t = 0 is given by

\hat{ρ} (0) = | g (⟩ ⟨) g | {\hat{ρ}}_{b},

(13)

where ${\hat{ρ}}_{b} = e^{- β Ĥ_{b}} / T r_{b} {e^{- β Ĥ_{b}}}$ with β = 1/k_BT. The radiation is turned on at t = 0. Thus, for t ≥ 0, the total Hamiltonian becomes Eq. (9). The corresponding time evolution operator for the material is

{\hat{U}}_{T} (t, 0) = \exp_{(+)} \{- \frac{i}{ℏ} \int_{0}^{t} d t^{'} Ĥ_{T} (t^{'})\},

(14)

and the density operator at time t is given by

\hat{ρ} (t) = {\hat{U}}_{T} (t, 0) \hat{ρ} (0) {\hat{U}}_{T}^{†} (t, 0) .

(15)

The linear absorption spectroscopy can then be expressed as the steady state limit of the time derivative of the lowest (the second) order approximation for the exciton population calculated from the above density operator. Appendix A provides detailed expressions of the corresponding perturbative expansion.

In the presence of time dependent fluctuations for E_g(t) and E_e,0(t) as assumed here, additional averaging over the fluctuations is necessary. Namely, the linear absorption line shape can be approximated as an average of Eq. (A4) over time dependent fluctuations, which is then divided by a normalization factor, ℏ²/(2π|A|²). Thus,

\begin{matrix} I_{a b s} (ω_{r}) & = & \frac{ℏ^{2}}{2 π | A |^{2}} lim_{t \to \infty} ⟨\frac{d}{d t} P_{e}^{(2)} (t)⟩ \\ = & \frac{1}{π} R e lim_{t \to \infty} \int_{0}^{t} d τ e^{i ω_{r} τ} ⟨) U_{g}^{†} (t, t - τ) U_{e, 0} (t, t - τ) ⟩) \\ \times {T r}_{e b} \{e^{i Ĥ_{b} τ / ℏ} Û_{e b} (τ) | D_{e_{r}} (⟩ ⟨) D_{e_{r}} | {\hat{ρ}}_{b}\}, \end{matrix}

(16)

where Tr_eb represents trace over both the exciton states and the bath degrees of freedom, and ⟨⋯⟩ denotes averaging over time dependent fluctuations. Under the assumption that these fluctuations are stationary, the result of averaging is independent of t and can be simplified to

\begin{matrix} I_{a b s} (ω_{r}) & = & \frac{1}{π} R e \int_{0}^{\infty} d t e^{i ω_{r} t} ⟨) U_{g}^{†} (t, 0) U_{e, 0} (t, 0) ⟩) \\ \times {T r}_{e b} \{e^{i Ĥ_{b} t / ℏ} Û_{e b} (t) | D_{e_{r}} (⟩ ⟨) D_{e_{r}} | {\hat{ρ}}_{b}\} . \end{matrix}

(17)

Let us denote the operator inside the trace operation in the above equation, which represents the coherence between the excited and the ground electronic states, as follows:

{\hat{ρ}}_{e g} (t) = e^{i Ĥ_{b} t / ℏ} Û_{e b} (t) | D_{e_{r}} (⟩ ⟨) D_{e_{r}} | {\hat{ρ}}_{b} .

(18)

Let us also denote the trace of the above operator over the bath degrees of freedom as follows:

{\hat{σ}}_{e g} (t) = {T r}_{b} \{{\hat{ρ}}_{e g} (t)\} .

(19)

Then, Eq. (17) can be expressed as

I_{a b s} (ω_{r}) \approx \frac{1}{π} R e \int_{0}^{\infty} d t e^{i ω_{r} t} ⟨) U_{g}^{†} (t, 0) U_{e, 0} (t, 0) ⟩) {T r}_{e} {{\hat{σ}}_{e g} (t)},

(20)

where Tr_e represents the trace over the exciton space only (excluding the ground electronic state and the bath).

A formally exact time evolution equation for ${\hat{σ}}_{e g} (t)$ can be derived employing the standard projection operator formalism as detailed in Appendix B. When approximated up to the fourth order, one can obtain the following time-convolution equation for ${\hat{σ}}_{e g} (t)$ by taking trace of Eq. (B13) over the bath:

\begin{matrix} \frac{d}{d t} {\hat{σ}}_{e g} (t) & \approx & - \frac{i}{ℏ} Ĥ_{e} {\hat{σ}}_{e g} (t) - \frac{1}{ℏ^{2}} \int_{0}^{t} d τ {\hat{C}}_{2} (t - τ) {\hat{σ}}_{e g} (τ) \\ + \frac{i}{ℏ^{3}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ {\hat{C}}_{3} (t - τ^{'}, τ^{'} - τ) {\hat{σ}}_{e g} (τ) \\ + \frac{1}{ℏ^{4}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ^{''} \int_{0}^{τ^{''}} d τ \\ \times {\hat{C}}_{4} (t - τ^{'}, τ^{'} - τ^{''}, τ^{''} - τ) {\hat{σ}}_{e g} (τ), \end{matrix}

(21)

where

{\hat{C}}_{2} (t) = {T r}_{b} \{Ĥ_{e b} (t) Û_{e} (t) Ĥ_{e b} {\hat{ρ}}_{b}\},

(22)

{\hat{C}}_{3} (t, t^{'}) = {T r}_{b} \{Ĥ_{e b} (t + t^{'}) Û_{e} (t) Ĥ_{e b} (t^{'}) Û_{e} (t^{'}) Ĥ_{e b} {\hat{ρ}}_{b}\},

(23)

\begin{matrix} {\hat{C}}_{4} (t, t^{'}, t^{''}) & = & {T r}_{b} \{Ĥ_{e b} (t + t^{'} + t^{''}) Û_{e} (t) Ĥ_{e b} (t^{'} + t^{''})) \\ \times Û_{e} (t^{'}) Ĥ_{e b} (t^{''}) Û_{e} (t^{''}) Ĥ_{e b} {\hat{ρ}}_{b} (\} \\ - {T r}_{b} \{Ĥ_{e b} (t) Û_{e} (t) Ĥ_{e b} {\hat{ρ}}_{b}\} Û_{e} (t^{'}) \\ \times {T r}_{b} \{Ĥ_{e b} (t^{''}) Û_{e} (t^{''}) Ĥ_{e b} {\hat{ρ}}_{b}\} . \end{matrix}

(24)

In the above expressions, Ĥ_eb(t) is the exciton-bath Hamiltonian represented in the partial interaction picture (with respect to the bath) representation and is given by

Ĥ_{e b} (t) = e^{i Ĥ_{b} t / ℏ} Ĥ_{e b} e^{- Ĥ_{b} t / ℏ} .

(25)

Let us introduce the eigenstates and eigenvalues of Ĥ_e as follows:

Ĥ_{e} | φ_{p} (⟩ = E_{p} | φ_{p} (⟩, p = 1, \dots, N_{c} .

(26)

Then, in terms of the unitary transformation $U$ with matrix elements, $U_{j p} = ⟨) s_{j} | φ_{p} (⟩$ ⁠, the site excitation state can be expressed as

| s_{j} (⟩ = \sum_{p = 1}^{N_{c}} U_{j p}^{*} | φ_{p} (⟩ .

(27)

When represented in the basis of |φ_p⟩’s, the transition dipole state |D_e⟩ defined by Eq. (11) can also be expressed as

| D_{e} (⟩ = \sum_{p} \sum_{j} U_{j p}^{*} (μ_{j} \cdot e) | φ_{p} (⟩ = \sum_{p = 1}^{N_{c}} ({\tilde{μ}}_{p} \cdot e) | φ_{p} (⟩,

(28)

where ${\tilde{μ}}_{p} = \sum_{j = 1}^{N_{c}} U_{j p}^{*} μ_{j}$ is the transition dipole for the excitation to |φ_p⟩ from the ground state |g⟩.

Employing Eq. (27) in Eq. (3), the exciton-bath Hamiltonian in the basis of exciton eigenstates can be expressed as

Ĥ_{e b} = \sum_{p = 1}^{N_{c}} \sum_{q = 1}^{N_{c}} \sum_{j = 1}^{N_{c}} \sum_{k = 1}^{N_{c}} U_{j p}^{*} {\hat{B}}_{j k} U_{k q} | φ_{p} (⟩ ⟨) φ_{q} | = \sum_{p = 1}^{N_{c}} \sum_{q = 1}^{N_{c}} {\hat{B}}_{p q} | φ_{p} (⟩ ⟨) φ_{q} |,

(29)

where the second equality serves as the definition of ${\hat{B}}_{p q}$ ⁠. Then, Eqs. (22)–(24) can be expressed as

{\hat{C}}_{2} (t) = \sum_{p} \sum_{q} | φ_{p} (⟩ ⟨) φ_{q} | \sum_{p_{1}} e^{- i E_{p_{1}} t / ℏ} {T r}_{b} \{{\hat{B}}_{p p_{1}} (t) {\hat{B}}_{p_{1} q} {\hat{ρ}}_{b}\},

(30)

\begin{matrix} {\hat{C}}_{3} (t, t^{'}) & = & \sum_{p} \sum_{q} | φ_{p} (⟩ ⟨) φ_{q} | \sum_{p_{1}} \sum_{p_{2}} e^{- i (E_{p_{1}} t + E_{p_{2}} t^{'}) / ℏ} \\ \times {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'}) {\hat{B}}_{p_{1} p_{2}} (t^{'}) {\hat{B}}_{p_{2} q} {\hat{ρ}}_{b}\}, \end{matrix}

(31)

\begin{matrix} {\hat{C}}_{4} (t, t^{'}, t^{''}) & = & \sum_{p} \sum_{q} | φ_{p} (⟩ ⟨) φ_{q} | \sum_{p_{1}} \sum_{p_{2}} \sum_{p_{3}} e^{- i (E_{p_{1}} t + E_{p_{2}} t^{'} + E_{p_{3}} t^{''}) / ℏ} \\ \times () {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'} + t^{''}) {\hat{B}}_{p_{1} p_{2}} (t^{'} + t^{''}) {\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} \\ - {T r}_{b} \{{\hat{B}}_{p p_{1}} (t) {\hat{B}}_{p_{1} p_{2}} {\hat{ρ}}_{b}\} {T r}_{b} \{{\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\})) . \end{matrix}

(32)

The time dependence of $⟨) U_{g}^{†} (t, 0) U_{e, 0} (t, 0) ⟩)$ in Eq. (17) can vary depending on the specific nature of external conditions and environments. In this work, let us consider the following simple case:

⟨) U_{g}^{†} (t, 0) U_{e, 0} (t, 0) ⟩) = e^{- i ⟨ E_{e g} ⟩ t / ℏ - Γ_{f l} t},

(33)

where Γ_fl is related to the second cumulant of the fluctuation of E_e,0(t) − E_g(t) in the fast modulation limit. Note that the source of Γ_fl is different from the spontaneous decay, but it has the same effect as the latter on the line shape expression. The assumption here is that Γ_fl is much larger than the spontaneous decay rate. Then, introducing

Ω = ω_{r} - ⟨ E_{e g} ⟩ / ℏ + i Γ_{f l},

(34)

the line shape can be expressed as follows:

I_{a b s} (ω_{r}) \approx \frac{1}{π} R e {T r}_{e} \{) \int_{0}^{\infty} d t e^{i Ω t} {\hat{σ}}_{e g} (t) \}) .

(35)

Multiplying Eq. (21) with e^iΩt and then integrating it over time t,

\begin{matrix} {(e^{i Ω t} {\hat{σ}}_{e g} (t)|}_{0}^{\infty} - i Ω \int_{0}^{\infty} d t e^{i Ω t} {\hat{σ}}_{e g} (t) \\ = (- \frac{i}{ℏ} Ĥ_{e x} - {\hat{K}}_{e b} (Ω)) \int_{0}^{\infty} d t e^{i Ω t} {\hat{σ}}_{e g} (t), \end{matrix}

(36)

where

{\hat{K}}_{e b} (Ω) \approx {\hat{K}}_{e b}^{(2)} (Ω) + {\hat{K}}_{e b}^{(3)} (Ω) + {\hat{K}}_{e b}^{(4)} (Ω),

(37)

with

{\hat{K}}_{e b}^{(2)} (Ω) = \frac{1}{ℏ^{2}} \int_{0}^{\infty} d t e^{i Ω t} {\hat{C}}_{2} (t),

(38)

\begin{matrix} {\hat{K}}_{e b}^{(3)} (Ω) = - \frac{i}{ℏ^{3}} \int_{0}^{\infty} d t \int_{0}^{\infty} d t^{'} e^{i Ω (t + t^{'})} {\hat{C}}_{3} (t, t^{'}), \\ {\hat{K}}_{e b}^{(4)} (Ω) = - \frac{1}{ℏ^{4}} \int_{0}^{\infty} d t \int_{0}^{\infty} d t^{'} \int_{0}^{\infty} d t^{''} e^{i Ω (t + t^{'} + t^{''})} {\hat{C}}_{4} (t, t^{'}, t^{''}) . \end{matrix}

(39)

Employing Eqs. (30)–(32), the matrix elements of the above operators in the basis of the eigenstates of Ĥ_e can be shown to be

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(2)} (Ω) | φ_{q} (⟩ & = & \frac{1}{ℏ^{2}} \sum_{p_{1}} \int_{0}^{\infty} d t e^{i (Ω - E_{p_{1}} / ℏ) t} \\ \times {T r}_{b} \{{\hat{B}}_{p p_{1}} (t) {\hat{B}}_{p_{1} q} {\hat{ρ}}_{b}\}, \end{matrix}

(40)

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(3)} (Ω) | φ_{q} (⟩ & = & - \frac{i}{ℏ^{3}} \sum_{p_{1}} \sum_{p_{2}} \int_{0}^{\infty} d t e^{i (Ω - E_{p_{1}} / ℏ) t} \\ \times \int_{0}^{\infty} d t^{'} e^{i (Ω - E_{p_{2}} / ℏ) t^{'}} \\ \times {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'}) {\hat{B}}_{p_{1} p_{2}} (t^{'}) {\hat{B}}_{p_{2} q} {\hat{ρ}}_{b}\}, \end{matrix}

(41)

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(4)} (Ω) | φ_{q} (⟩ & = & - \frac{1}{ℏ^{4}} \sum_{p_{1}} \sum_{p_{2}} \sum_{p_{3}} \int_{0}^{\infty} d t e^{i (Ω - E_{p_{1}} / ℏ) t} \\ \times \int_{0}^{\infty} d t^{'} e^{i (Ω - E_{p_{2}} / ℏ) t} \int_{0}^{\infty} d t^{''} e^{i (Ω - E_{p_{3}} / ℏ) t^{''}} \\ \times ({T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'} + t^{''}) {\hat{B}}_{p_{1} p_{2}} (t^{'} + t^{''}))) \\ ((\times {\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} - {T r}_{b} \{{\hat{B}}_{p p_{1}} (t) {\hat{B}}_{p_{1} p_{2}} {\hat{ρ}}_{b}\}) \\ (\times {T r}_{b} \{{\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\}) . \end{matrix}

(42)

In Eq. (36), ${lim}_{t \to \infty} e^{i Ω t} {\hat{σ}}_{e g} (t) = 0$ because of the positive value of Γ_fl. Thus, it can be solved as follows:

\int_{0}^{\infty} d t e^{i Ω t} {\hat{σ}}_{e g} (t) = i {[Ω - Ĥ_{e} / ℏ + i {\hat{K}}_{e b} (Ω)]}^{- 1} | D_{e} (⟩ ⟨) D_{e} | .

(43)

Then, the absorption line shape defined by Eq. (35) can be expressed as

I_{a b s} (Ω) = - \frac{1}{π} I m ⟨) D_{e} | {[Ω - Ĥ_{e} / ℏ + i {\hat{K}}_{e b} (Ω)]}^{- 1} | D_{e} (⟩,

(44)

where the fact that ${\hat{1}}_{e} = \sum_{j} | s_{j} (⟩ ⟨) s_{j} |$ forms an effective identity operator in the exciton space has been used.

While direct numerical evaluation of Eq. (44) through inversion of the complex matrix is not difficult, it is instructive to consider its secular approximation, which amounts to assuming that ${\hat{K}}_{e b} (Ω)$ is diagonal in the basis of the eigenstates of Ĥ_e as follows:

{\hat{K}}_{e b} (Ω) \approx \sum_{p} (K_{e b, p}^{R} (Ω) + i K_{e b, p}^{I} (Ω)) | φ_{p} (⟩ ⟨) φ_{p} |,

(45)

where

K_{e b, p}^{R} (Ω) = R e ⟨) φ_{p} | {\hat{K}}_{e b} (Ω) | φ_{p} (⟩,

(46)

K_{e b, p}^{I} (Ω) = I m ⟨) φ_{p} | {\hat{K}}_{e b} (Ω) | φ_{p} (⟩ .

(47)

Then, the line shape expression simplifies to

I_{a b s} (ω_{r}) \approx \frac{1}{π} \sum_{p = 1}^{N_{c}} \frac{| {\tilde{μ}}_{p} \cdot e_{r} |^{2} K_{e b, p}^{R} (Ω)}{{[Ω - K_{e b, p}^{I} (Ω)]}^{2} + K_{e b, p}^{R} {(Ω)}^{2}} .

(48)

In the above expression, $K_{e b, p}^{I} (Ω)$ represents the shift of the excitation energy due to interaction with the bath and $K_{e b, p}^{R} (Ω)$ represents the broadening due to bath relaxation. It is important to note that the effects of time dependent fluctuations are incorporated into both terms through the imaginary component, Γ_fl, of Ω.

III. LINE SHAPE EXPRESSIONS FOR A LINEARLY COUPLED HARMONIC OSCILLATOR BATH MODEL

Consider the well-known bath model of linearly coupled harmonic oscillators with the following Hamiltonians:

Ĥ_{e b} = \sum_{j} \sum_{n} ℏ ω_{n} g_{n j} ({\hat{b}}_{n} + {\hat{b}}_{n}^{†}) | s_{j} (⟩ ⟨) s_{j} |,

(49)

Ĥ_{b} = \sum_{n} ℏ ω_{n} ({\hat{b}}_{n}^{†} {\hat{b}}_{n} + \frac{1}{2}),

(50)

where ${\hat{b}}_{n}^{†}$ and ${\hat{b}}_{n}$ are creation and annihilation operators of each bath mode. For this case, the bath operator defined by Eq. (29) becomes

{\hat{B}}_{p q} = \sum_{j} \sum_{n} ℏ ω_{n} g_{n j} U_{j p}^{*} U_{j q} ({\hat{b}}_{n} + {\hat{b}}_{n}^{†}) .

(51)

Then, Eq. (40) for the present model is expressed as

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(2)} (Ω) | φ_{q} (⟩ \\ = \sum_{p_{1}} \sum_{j} \sum_{j^{'}} \sum_{n} ω_{n}^{2} g_{n j} g_{n j^{'}} U_{j p}^{*} U_{j p_{1}} U_{j^{'} p_{1}}^{*} U_{j^{'} q} \\ \times \int_{0}^{\infty} d t e^{i (Ω - E_{p_{1}} / ℏ) t} () ⟨ N_{n} + 1 ⟩ e^{- i ω_{n} t} + ⟨ N_{n} ⟩ e^{i ω_{n} t})), \end{matrix}

(52)

where $⟨ N_{n} ⟩ = 1 / (e^{β ℏ ω_{n}} - 1)$ ⁠, which is the thermal average of the vibrational quanta of mode n. Introducing the following spectral density,

J_{j j^{'}} (ω) = π ℏ \sum_{n} δ (ω - ω_{n}) ω_{n}^{2} g_{j n} g_{j n^{'}},

(53)

and doing the integration over t, Eq. (52) can be expressed as

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(2)} (Ω) | φ_{q} (⟩ & = & \frac{i}{π ℏ} \sum_{j} \sum_{j^{'}} \sum_{p_{1}} U_{j p}^{*} U_{j p_{1}} U_{j^{'} p_{1}}^{*} U_{j^{'} q} \\ \times \int_{0}^{\infty} d ω J_{j j^{'}} (ω) \{) \frac{1}{(1 - e^{- β ℏ ω})} \frac{1}{Ω - E_{p_{1}} / ℏ - ω} \\ + \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{1}{(Ω - E_{p_{1}} / ℏ + ω)} \}) . \end{matrix}

(54)

Let us define

Ω_{p, r} = ω_{r} - \frac{⟨ E_{e g} ⟩ + E_{p}}{ℏ} .

(55)

Then, $Ω - E_{p} / ℏ = Ω_{p, r} + i Γ_{f l}$ ⁠. Employing this expression, Eq. (54) can be expressed as

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(2)} (Ω) | φ_{q} (⟩ \\ = \frac{1}{π ℏ} \sum_{j} \sum_{j^{'}} \sum_{p_{1}} U_{j p}^{*} U_{j p_{1}} U_{j^{'} p_{1}}^{*} U_{j^{'} q} \int_{0}^{\infty} d ω J_{j j^{'}} (ω) \{) \frac{1}{(1 - e^{- β ℏ ω})} \\ \times \frac{Γ_{f l} + i (Ω_{p_{1}, r} - ω)}{{(Ω_{p_{1}, r} - ω)}^{2} + Γ_{f l}^{2}} + \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{Γ_{f l} + i (Ω_{p_{1}, r} + ω)}{{(Ω_{p_{1}, r} + ω)}^{2} + Γ_{f l}^{2}} \}) . \end{matrix}

(56)

For the present model, the third order bath correlation term $⟨) φ_{p} | {\hat{K}}_{e b}^{(3)} (Ω) | φ_{q} (⟩ = 0$ ⁠.

The fourth order bath correlation function in Eq. (42) for the present model can be simplified employing Wick’s theorem, which can also be proven explicitly as shown in Appendix C. Employing Eq. (C9) and conducting integrations over time, the following expression can be obtained:

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(4)} (Ω) | φ_{q} (⟩ & = & \frac{i}{{(π ℏ)}^{2}} \sum_{p_{1}, p_{2}, p_{3}} \sum_{j_{1}, j_{2}, j_{3}, j_{4}} U_{j_{1} p}^{*} U_{j_{1} p_{1}} U_{j_{2} p_{1}}^{*} U_{j_{2} p_{2}} U_{j_{3} p_{2}}^{*} U_{j_{3} p_{3}} U_{j_{4} p_{3}}^{*} U_{j_{4} q} \\ \times [) \int_{0}^{\infty} d ω \frac{J_{j_{1} j_{3}} (ω)}{(1 - e^{- β ℏ ω})} \int_{0}^{\infty} d ω^{'} \frac{J_{j_{2} j_{4}} (ω^{'})}{(1 - e^{- β ℏ ω^{'}})} \{\frac{1}{(Ω_{p 1} - ω) (Ω_{p 2} - ω - ω^{'}) (Ω_{p 3} - ω^{'})}) \\ (+ \frac{e^{- β ℏ ω - β ℏ ω^{'}}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω + ω^{'}) (Ω_{p 3} + ω^{'})} + \frac{e^{- β ℏ ω^{'}}}{(Ω_{p 1} - ω) (Ω_{p 2} - ω + ω^{'}) (Ω_{p 3} + ω^{'})} + \frac{e^{- β ℏ ω}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω - ω^{'}) (Ω_{p 3} - ω^{'})}\} \\ + \int_{0}^{\infty} d ω \frac{J_{j_{1} j_{4}} (ω)}{(1 - e^{- β ℏ ω})} \int_{0}^{\infty} d ω^{'} \frac{J_{j_{2} j_{3}} (ω^{'})}{(1 - e^{- β ℏ ω^{'}})} \{\frac{1}{(Ω_{p 1} - ω) (Ω_{p 2} - ω - ω^{'}) (Ω_{p 3} - ω)} + \frac{e^{- β ℏ ω - β ℏ ω^{'}}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω + ω^{'}) (Ω_{p 3} + ω)}) \\ (+ \frac{e^{- β ℏ ω^{'}}}{(Ω_{p 1} - ω) (Ω_{p 2} - ω + ω^{'}) (Ω_{p 3} - ω)} + \frac{e^{- β ℏ ω}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω - ω^{'}) (Ω_{p 3} + ω)}\}]) . \end{matrix}

(57)

In the above expression, Ω_p = Ω_p,r + iΓ_fl with Ω_p,r as defined by Eq. (55).

IV. MODEL CALCULATIONS

For the numerical test, we consider a linear chain of identical chromophores with only nearest neighbor interactions. In the absence of any disorder or dynamic fluctuations, this model has exact solutions for the single exciton states and energies. It is assumed that the excited states in each chromophore unit are well separated so that we can focus on only one excited state with energy E. The nearest neighbor electronic coupling is denoted as J. Then, the Hamiltonian representing the single exciton states that are formed by all the states with energy E, i.e., |s_j⟩ with j = 1, …, N_c, is given by

Ĥ_{e} = E \sum_{j = 1}^{N_{c}} | s_{j} (⟩ ⟨) s_{j} | + J \sum_{j = 1}^{N_{c} - 1} (| s_{j} (⟩ ⟨) s_{j + 1} | + | s_{j + 1} (⟩ ⟨) s_{j} |) .

(58)

The eigenstates and eigenvalues of this Hamiltonian are as follows:

E_{p} = E + 2 J \cos (\frac{p π}{N_{c} + 1}),

(59)

| φ_{p} (⟩ = {(\frac{2}{N_{c} + 1})}^{1 / 2} \sum_{j = 1}^{N_{c}} \sin (\frac{p j π}{N_{c} + 1}) | s_{j} (⟩,

(60)

where p = 1, …, N_c. Assume that the jth molecular unit has the following transition dipole moment:

μ_{j} = μ (\sin θ \cos (ϕ_{0} + j δ ϕ), \sin θ \sin (ϕ_{0} + j δ ϕ), \cos θ),

(61)

where the z-direction has been defined as the direction of the linear chain. Thus, θ is the axial angle and ϕ₀ + jδϕ is the azimuthal angle of the transition dipole μ_j. For this, the components of the transition dipole vector for each eigenstate can be shown to be

\begin{matrix} {\tilde{μ}}_{p, x} & = & \frac{\sin θ}{{[8 (N_{c} + 1)]}^{1 / 2}} \{\frac{\cos (ϕ_{p}^{+} + ϕ_{0})}{\sin (ϕ_{p}^{+})} + \frac{\cos (ϕ_{p}^{-} - ϕ_{0})}{\sin (ϕ_{p}^{-})}) \\ - {(- 1)}^{p} \frac{\cos (ϕ_{p}^{-} - N_{c} δ ϕ - ϕ_{0})}{\sin (ϕ_{p}^{+})} \\ (- {(- 1)}^{p} \frac{\cos (ϕ_{p}^{+} + N δ ϕ + ϕ_{0})}{\sin (ϕ_{p}^{-})}\}, \end{matrix}

(62)

\begin{matrix} {\tilde{μ}}_{p, y} & = & \frac{\sin θ}{{[8 (N_{c} + 1)]}^{1 / 2}} \{\frac{\sin (ϕ_{p}^{+} + ϕ_{0})}{\sin (ϕ_{p}^{+})} - \frac{\sin (ϕ_{p}^{-} - ϕ_{0})}{\sin (ϕ_{p}^{-})}) \\ (+ {(- 1)}^{p} \frac{\cos (ϕ_{p}^{-} - N_{c} δ ϕ - ϕ_{0})}{\sin (ϕ_{p}^{+})} - {(- 1)}^{p} \frac{\sin (ϕ_{p}^{+} + N_{c} δ ϕ + ϕ_{0})}{\sin (ϕ_{p}^{-})}\}, \end{matrix}

(63)

{\tilde{μ}}_{p, z} = \frac{(1 - {(- 1)}^{p}) \cos θ}{{[2 (N_{c} + 1)]}^{1 / 2}} \cot (\frac{p π}{2 (N_{c} + 1)}),

(64)

where

ϕ_{p}^{+} = \frac{1}{2} (\frac{p π}{(N_{c} + 1)} + δ ϕ),

(65)

ϕ_{p}^{-} = \frac{1}{2} (\frac{p π}{(N_{c} + 1)} - δ ϕ) .

(66)

It is assumed that the site excitation has the same local modes and there is no mode commonly coupled to different site excitations. Thus,

J_{j j^{'}} (ω) = δ_{j j^{'}} J (ω) .

(67)

Then, introduce the following coefficient:

\begin{matrix} W_{p p^{'} p^{''} q} & = & \sum_{j} U_{j p}^{*} U_{j p^{'}} U_{j p^{''}}^{*} U_{j q} \\ = & \frac{4}{{(N_{c} + 1)}^{2}} \sum_{j = 1}^{N_{c}} \sin (\frac{p j π}{N_{c} + 1}) \sin (\frac{p^{'} j π}{N_{c} + 1}) \\ \times \sin (\frac{p^{''} j π}{N_{c} + 1}) \sin (\frac{q j π}{N_{c} + 1}), \end{matrix}

(68)

which becomes the participation ratio for $p = p^{'} = p^{''} = q$ ⁠, Eq. (56) can then be expressed as

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(2)} (Ω) | φ_{q} (⟩ & = & \frac{1}{π ℏ} \sum_{p_{1}} W_{p p_{1} p_{1} q} \int_{0}^{\infty} d ω J (ω) \\ \times \{) \frac{1}{(1 - e^{- β ℏ ω})} \frac{Γ_{f l} + i (Ω_{p_{1}, r} - ω)}{{(Ω_{p_{1}, r} - ω)}^{2} + Γ_{f l}^{2}} \\ + \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{Γ_{f l} + i (Ω_{p_{1}, r} + ω)}{{(Ω_{p_{1}, r} + ω)}^{2} + Γ_{f l}^{2}} \}) . \end{matrix}

(69)

Similarly, Eq. (57) simplifies to

\begin{matrix} ⟨) φ_{p} | {\hat{K}}_{e b}^{(4)} (Ω) | φ_{q} (⟩ & = & \frac{i}{{(π ℏ)}^{2}} \int_{0}^{\infty} d ω \frac{J (ω)}{(1 - e^{- β ℏ ω})} \int_{0}^{\infty} d ω^{'} \frac{J (ω^{'})}{(1 - e^{- β ℏ ω^{'}})} \sum_{p_{1}, p_{2}, p_{3}} [W_{p p_{1} p_{2} p_{3}} W_{p_{1} p_{2} p_{3} q} \{\frac{1}{(Ω_{p 1} - ω) (Ω_{p 2} - ω - ω^{'}) (Ω_{p 3} - ω^{'})})) \\ ((+ \frac{e^{- β ℏ ω - β ℏ ω^{'}}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω + ω^{'}) (Ω_{p 3} + ω^{'})} + \frac{e^{- β ℏ ω^{'}}}{(Ω_{p 1} - ω) (Ω_{p 2} - ω + ω^{'}) (Ω_{p 3} + ω^{'})} + \frac{e^{- β ℏ ω}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω - ω^{'}) (Ω_{p 3} - ω^{'})})\} \\ (+ W_{p p_{1} p_{3} q} W_{p_{1} p_{2} p_{2} p_{3}} \{\frac{1}{(Ω_{p 1} - ω) (Ω_{p 2} - ω - ω^{'}) (Ω_{p 3} - ω)})) + \frac{e^{- β ℏ ω - β ℏ ω^{'}}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω + ω^{'}) (Ω_{p 3} + ω)} \\ (+ \frac{e^{- β ℏ ω^{'}}}{(Ω_{p 1} - ω) (Ω_{p 2} - ω + ω^{'}) (Ω_{p 3} - ω)} (+ \frac{e^{- β ℏ ω}}{(Ω_{p 1} + ω) (Ω_{p 2} + ω - ω^{'}) (Ω_{p 3} + ω)}\}] . \end{matrix}

(70)

Assuming that J can be modeled by a sum of an effective transition dipole interaction and an additional term J₀ representing the remaining sum of exchange and higher order Coulomb interactions, the nearest neighbor electronic coupling can be expressed as

J = J_{0} + \frac{μ^{2}}{ϵ d^{3}} (\sin^{2} θ \cos δ ϕ - 2 \cos^{2} θ) .

(71)

All the parameters introduced below are defined in the Lorentz-Heaviside units of ℏ = c = 1. In addition, Å is defined as the unit length and is assumed to be equal to d, the nearest neighbor chromophore distance. In these units, the value of the transition dipole moment μ = 10, which corresponds to about 4.09 D, will be used for all the calculations. For convenience, it is assumed that ϵ = 1.

While various values of the angles θ, δϕ, and ϕ₀ can be employed, only the case with θ = 60°, δϕ = 24°, and ϕ₀ = 0 is considered here. For this choice of θ and δϕ, the nearest neighbor transition dipole coupling is 18.5, and the second nearest transition dipole interaction almost vanishes, which justifies the assumption of only nearest neighbor electronic couplings. The bath spectral density is assumed to be a super-Ohmic spectral form as follows: $J (ω) = π η (ω^{3} / ω_{c}^{2}) e^{- ω / ω_{c}}$ with ω_c = 100.

Line shapes based on both the 2nd and 4th order approximations were calculated for two cases of J₀ = −200 and 200 with N_c = 3 and 6 and assuming Γ_fl = 10. Two different values of η = 0.2 and η = 0.5 were considered at k_BT = 200, which is close to room temperature. Also compared are the line shapes calculated by the following diagonal approximation in the exciton basis:

\begin{matrix} I_{a b s}^{d} & = & \frac{1}{π} R e \{) \sum_{p} {|e \cdot {\tilde{μ}}_{p}|}^{2} \int_{0}^{\infty} d t e^{i t (Ω_{p} + I_{p} λ / ℏ) t} \\ \times \exp \{) - I_{p} G_{r} (t) - i I_{p} G_{i} (t) \}) \}), \end{matrix}

(72)

where I_p = W_pppp, the participation ratio of the exciton state p, and

\begin{aligned} G_{r} (t) & = \frac{1}{π ℏ} \int_{0}^{\infty} d ω \frac{J (ω)}{ω^{2}} \coth (\frac{β ℏ ω}{2}) (1 - \cos (ω t)), \\ G_{i} (t) & = \frac{1}{π ℏ} \int_{0}^{\infty} d ω \frac{J (ω)}{ω^{2}} \sin (ω t), \end{aligned}

(73)

λ = \frac{1}{π} \int_{0}^{\infty} d ω \frac{J (ω)}{ω} .

(74)

Within this approximation, the diagonal components of the exciton-bath are included up to the infinite order. This is simpler than the modified Redfield equation approach^25,27 but serves as a good reference for a clear physical understanding of the features of QME-based line shapes. The resulting line shapes averaged over all polarization directions are provided in Figs. 1 and 2.

FIG. 1.

View large Download slide

Absorption line shapes for J₀ = −200 and Γ_fl = 10 at k_BT = 200 in the units of ℏ = c = 1, and Å as unit length, for η = 0.2 and 0.5 and for N_c = 3 and 6. QME2 represents the 2nd order QME calculation result, and QME4 represents the 4th order QME calculation result. “Diagonal approx.” represents the approximation, Eq. (72), which includes only the diagonal components of the exciton-bath coupling in the exciton basis.

FIG. 2.

View large Download slide

Absorption line shapes for J₀ = 200 and k_BT = 200 in the units of ℏ = c = 1 and Å as unit length. Other conventions are the same as Fig. 1.

The case with J₀ = −200 shows J-aggregate-like behavior, whereas J₀ = 200 has more H-aggregate-like feature. Thus, in the former case, lower exciton states have dominant oscillator strengths, whereas the opposite situation occurs for the latter case. The peak positions for all of the 2nd and 4th order QMEs and the diagonal exciton-bath coupling approximation are fairly close for N_c = 3 and η = 0.2. However, the linewidths for higher exciton states based on the QME methods are significantly broader than those for the diagonal exciton-bath approximation. This is because the former can account for transitions between different excitons through vibronic couplings. For η = 0.5, there are significant discrepancies in both peak positions and linewidths.

As is clear from Figs. 1 and 2, the 4th order results add additional features to the 2nd ones that amount to higher order exciton-phonon effects. However, it is also apparent that consideration of up to 4th order terms is not sufficient for developing a reliable method for a broad range of parameters because of the breakdown of the positive definiteness. While the line shapes based on the fourth order calculation remain positive for η = 0.5, a test for a larger value of η showed that the line shape can become negative at some values of the frequency.

Figures 3 and 4 compare the full 4th order QME line shape expressions with the secular approximation [Eq. (48)]. For η = 0.5, the two are very close. For the case η = 0.2, which is not shown here, the two are in fact almost identical. However, for η = 1, a significant difference can be seen between the two. More importantly, there are narrow ranges of frequency where the line shape becomes negative for both full 4th order and the secular approximation. A simple remedy for this is to use the following extrapolation within the secular approximation:

K_{e b, p}^{R} (Ω) \approx K_{e b, p}^{(2), R} (Ω) \exp \{\frac{K_{e b, p}^{(4), R} (Ω)}{K_{e b, p}^{(2), R} (Ω)}\},

(75)

where the superscript R represents the real part. Figure 3 shows that this extrapolation results in positive line shapes. This shows the importance of a closure or extrapolation scheme in using perturbative approximations, for which Padé^40,47,51 resummation and Landau-Zenner⁴⁹ resummation serve as good candidates.

FIG. 3.

View large Download slide

Comparison of full 4th order QME line shapes for J₀ = −200 with the secular approximation, Eq. (48), and also with the latter with an additional extrapolation, Eq. (75). Both cases of N_c = 3 and N_c = 6 for η = 0.5 and 1 are shown. All other parameters are the same as in Fig. 1.

FIG. 4.

View large Download slide

Comparison of full 4th order QME line shapes for J₀ = 200 with the secular approximation, Eq. (48), and also with the latter with an additional extrapolation. Equation (75). Both cases of N_c = 3 and N_c = 6 for η = 0.5 and 1 are shown. All other parameters are the same as for Fig. 2.

V. CONCLUSION

This work has provided general and detailed 4th order expressions of the absorption line shape for generic exciton-bath models by going through explicit calculation of all the fourth order terms within the time-convolution QME. The resulting expressions require additional integration in the frequency domain but are otherwise easy to calculate. The expressions also do not make any assumptions regarding the forms of the spectral density.

Results of calculation for a simple model of linear chain of chromophores demonstrate the feasibility of including the fourth order terms in general and provide insights into the features they can represent in the line shape expression. While it is clear that the 4th order terms capture new multiphonon terms missing in the 2nd order expressions, it is also seen that the inclusion of the higher order terms tends to make line shapes more prone to be nonpositive, which is unphysical, as the exciton-bath coupling becomes stronger. Therefore, for practical utility, an additional correction or extrapolation scheme that ensures positivity is necessary. Within the secular approximation, it is possible to find a simple remedy for this issue as can be seen from Figs. 3 and 4. However, a more careful numerical study is necessary in order to find a reliable extrapolation scheme that also results in correct behavior in the limit of strong exciton-bath coupling.

Numerical results presented here are as yet preliminary and do not serve as a satisfactory assessment or benchmark of the fourth order expressions that can address important pending issues. For example, it is already well established that the 2nd order QME expression is less satisfactory^26,27,34 than others for the case of Ohmic bath spectral density. However, the performance of higher order approximations and for other types of spectral densities⁵² remains unclear, which needs further numerical tests. The expressions presented in this work can serve as an important starting point for examining these issues. For this, test calculations for a broader range of models and spectral densities are needed. Benchmarking these results against numerically exact methods and other approximations^24,25,28 will be an important subject of future study.

ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (Grant No. CHE-1362926) during its initial stage of formal development and was completed with the support from the Office of Basic Energy Sciences, Department of Energy (Grant No. DE-SC0001393).

APPENDIX A: SECOND ORDER APPROXIMATION FOR THE EXCITON POPULATION

The lowest order approximation for the probability to find the material in the excited state can be obtained by expanding the time evolution operator up to the first order of Ĥ_er(t) as follows:

{\hat{U}}_{T}^{(1)} (t, 0) \approx Û_{m} (t, 0) - \frac{i}{ℏ} \int_{0}^{t} d t^{'} Û_{m} (t, t^{'}) Ĥ_{e r} (t^{'}) Û_{m} (t^{'}) .

(A1)

When approximated up to the second order of Ĥ_er(t), the probability to find the material in excited states at time t is

\begin{matrix} P_{e}^{(2)} (t) & = & \frac{1}{ℏ^{2}} \int_{0}^{t} d t^{'} \int_{0}^{t} d t^{''} \sum_{j = 1}^{N_{c}} {T r}_{b} \{⟨) s_{j} | U_{e, 0} (t, t^{'}) Û_{e b} (t, t^{'})) \\ \times Ĥ_{e r} (t^{'}) U_{g} (t^{'}, t^{''}) | g (⟩ ⟨) g | {\hat{ρ}}_{b} Ĥ_{e r} (t^{''}) \\ (\times Û_{e b}^{†} (t, t^{''}) U_{e, 0}^{*} (t, t^{''}) | s_{j} (⟩\}, \end{matrix}

(A2)

where $Û_{g} (t^{'}, t^{''})$ ⁠, U_e,0(t, t′), and Û_eb(t, t′) have been defined by Eqs. (5)–(7), respectively. Employing Eq. (10) for Ĥ_er(t) in the above expression, Eq. (A2) can be expressed as

\begin{matrix} P_{e}^{(2)} (t) & = & \frac{| A |^{2}}{ℏ^{2}} \int_{0}^{t} d t^{'} \int_{0}^{t} d t^{''} e^{- i ω_{r} (t^{'} - t^{''})} U_{g} (t^{'}, t^{''}) U_{e, 0} (t^{''}, t^{'}) \\ \times \sum_{j = 1}^{N_{c}} \sum_{j^{'} = 1}^{N_{c}} \sum_{j^{''} = 1}^{N_{c}} (μ_{j^{'}} \cdot e_{r}) (μ_{j^{''}}^{*} \cdot e_{r}^{*}) \\ \times {T r}_{b} \{⟨) s_{j} | Û_{e b} (t, t^{'}) | s_{j^{'}} (⟩ {\hat{ρ}}_{b} e^{i Ĥ_{b} (t^{''} - t^{'}) / ℏ} ⟨) s_{j^{''}} | Û_{e b}^{†} (t, t^{''}) | s_{j} (⟩\} . \end{matrix}

(A3)

Employing the definition of Eq. (11), Eq. (A3) can then be expressed as

\begin{matrix} P_{e}^{(2)} (t) & = & \frac{| A |^{2}}{ℏ^{2}} \int_{0}^{t} d t^{'} \int_{0}^{t} d t^{''} e^{- i ω_{r} (t^{'} - t^{''})} U_{g} (t^{'}, t^{''}) U_{e, 0} (t^{''}, t^{'}) \\ \times \sum_{j = 1}^{N_{c}} {T r}_{e b} \{e^{i Ĥ_{b} (t^{''} - t^{'}) / ℏ} Û_{e b} (t^{''}, t^{'}) | D_{e_{r}} (⟩ ⟨) D_{e_{r}} | {\hat{ρ}}_{b}\}, \end{matrix}

(A4)

where Tr_eb represents the trace over all the excitation states, |s_j⟩’s, and the bath. In the above expression, the invariance of the trace operation Tr_eb with respect to cyclic permutation of the operators, which are all defined in the subspace formed by the excited states and the bath, has been used.

APPENDIX B: EXACT FORMAL EXPRESSION AND FOURTH ORDER APPROXIMATION FOR THE TIME EVOLUTION OF PROJECTED DENSITY OPERATOR

Taking the time derivative of Eq. (18), it is straightforward to show that

\frac{d}{d t} {\hat{ρ}}_{e g} (t) = - \frac{i}{ℏ} (Ĥ_{e} + Ĥ_{e b} (t)) {\hat{ρ}}_{e g} (t),

(B1)

where Ĥ_eb(t) has been defined by Eq. (25). Let us introduce

Ĝ (t) = Û_{e}^{†} (t) {\hat{ρ}}_{e g} (t),

(B2)

where Û_e(t) has been defined by Eq. (8). Then,

\frac{d}{d t} Ĝ (t) = - \frac{i}{ℏ} Ĥ_{e b, I} (t) Ĝ (t),

(B3)

where

Ĥ_{e b, I} (t) = Û_{e}^{†} (t) H_{e b} (t) Û_{e} (t) .

(B4)

At this point, it is useful to introduce the most commonly used projection superoperator $P$ defined as

P (\cdot) = ρ_{b} {T r}_{b} \{(\cdot)\},

(B5)

where (·) represents the entire product of operators on the right-hand side of the projection superoperator. This satisfies the required idempotency property, $P^{2} = 1$ ⁠, and the following identity:

P Ĥ_{e} = Ĥ_{e} P .

(B6)

In addition, let us assume that the following condition is satisfied:

P Ĥ_{e b, I} (t) P = 0,

(B7)

which can always be satisfied by a proper definition of Ĥ_e(t). Then, by applying $P$ and $Q = 1 - P$ to Eq. (B3),

\frac{d}{d t} P Ĝ (t) = - \frac{i}{ℏ} P Ĥ_{e b, I} (t) Q Ĝ (t),

(B8)

\frac{d}{d t} Q Ĝ (t) = - \frac{i}{ℏ} Q Ĥ_{e b, I} (t) (Q + P) Ĝ (t) .

(B9)

The formal solution for $Q Ĝ (t)$ is as follows:

\begin{matrix} Q Ĝ (t) & = & \exp_{(+)} [- \frac{i}{ℏ} \int_{0}^{t} d τ Q Ĥ_{e b, I} (τ)] Q Ĝ (0) \\ - \frac{i}{ℏ} \int_{0}^{t} d τ \exp_{(+)} [- \frac{i}{ℏ} \int_{τ}^{t} d τ^{'} Q Ĥ_{e b, I} (τ^{'})] P Ĝ (τ) . \end{matrix}

(B10)

Inserting this into Eq. (B8), we obtain

\begin{matrix} \frac{d}{d t} P Ĝ (t) & = & - \frac{i}{ℏ} P Ĥ_{e b, I} (t) \exp_{(+)} [- \frac{i}{ℏ} \int_{0}^{t} d τ Q Ĥ_{e b, I} (τ)] Q Ĝ (0) \\ - \frac{1}{ℏ^{2}} \int_{0}^{t} d τ P Ĥ_{e b, I} (t) \exp_{(+)} [- \frac{i}{ℏ} \int_{τ}^{t} d τ^{'} Q Ĥ_{e b, I} (τ^{'})] \\ \times Q Ĥ_{e b, I} (τ) P Ĝ (τ) . \end{matrix}

(B11)

Combining the above equation with Eq. (B2) and using the property of Eq. (B6) lead to

\begin{matrix} \frac{d}{d t} P {\hat{ρ}}_{e g} (t) & = & - \frac{i}{ℏ} Ĥ_{e} P {\hat{ρ}}_{e g} (t) - \frac{1}{ℏ^{2}} Û_{e} (t) P Ĥ_{e b, I} (t) \int_{0}^{t} d τ \\ \times \exp_{(+)} [- \frac{i}{ℏ} \int_{τ}^{t} d τ^{'} Q Ĥ_{e b, I} (τ^{'})] Q Ĥ_{e b, I} (τ) Û_{e}^{†} (τ) P {\hat{ρ}}_{e g} (τ), \end{matrix}

(B12)

where the fact that $Q {\hat{ρ}}_{e g} (0) = 0$ has also been used. Up to the fourth order of Ĥ_eb,I(t), the above expression can be approximated as

\begin{matrix} \frac{d}{d t} P {\hat{ρ}}_{e g} (t) & \approx & - \frac{i}{ℏ} Ĥ_{e} P {\hat{ρ}}_{e g} (t) - \frac{1}{ℏ^{2}} P Ĥ_{e b} (t) [\int_{0}^{t} d τ Û_{e} (t - τ) Ĥ_{e b} (τ)) \\ - \frac{i}{ℏ} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ Û_{e} (t - τ^{'}) Ĥ_{e b} (τ^{'}) Û_{e} (τ^{'} - τ) Ĥ_{e b} (τ) \\ - \frac{1}{ℏ^{2}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ^{''} \int_{0}^{τ^{''}} d τ Û_{e} (t - τ^{'}) Ĥ_{e b} (τ^{'}) \\ \times (Û_{e} (τ^{'} - τ^{''}) Q Ĥ_{e b} (τ^{''}) Û_{e} (τ^{''} - τ) Ĥ_{e b} (τ)] P {\hat{ρ}}_{e g} (τ) . \end{matrix}

(B13)

The above equation is time nonlocal. A time local equation of the fourth order can also be derived by employing the following approximation:

P Ĝ (τ) \approx P Ĝ (t) + \frac{1}{ℏ^{2}} \int_{τ}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ^{''} P Ĥ_{e b, I} (τ^{'}) Q P Ĥ_{e b, I} (τ^{''}) P Ĝ (t) .

(B14)

The above relation can be derived by integrating the second order approximation for Eq. (B11) from τ to t and replacing $P Ĝ (τ^{''})$ within the integrand with $P Ĝ (t)$ ⁠, which does not affect the result up to the second order. Inserting the above approximation into the fourth order approximation for Eq. (B12), keeping only terms up to the fourth order of Ĥ_eb, and employing Eq. (B2), one can obtain the following expression:

\begin{matrix} \frac{d}{d t} P {\hat{ρ}}_{e g} (t) & \approx & - \frac{i}{ℏ} Ĥ_{e} P {\hat{ρ}}_{e g} (t) - \frac{1}{ℏ^{2}} [) \int_{0}^{t} d τ P Ĥ_{e b} (t) Û_{e} (t - τ) Ĥ_{e b} (τ) \\ - \frac{i}{ℏ} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ P Ĥ_{e b} (t) Û_{e} (t - τ^{'}) Ĥ_{e b} (τ^{'}) Û_{e} (τ^{'} - τ) Ĥ_{e b} (τ) \\ - \frac{1}{ℏ^{2}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ^{''} \int_{0}^{τ^{''}} d τ \{) P Ĥ_{e b} (t) Û_{e} (t - τ^{'}) Ĥ_{e b} (τ^{'}) Û_{e} (τ^{'} - τ^{''}) Q Ĥ_{e b} (τ^{''}) Û_{e} (τ^{''} - τ) Ĥ_{e b} (τ) \\ - P Ĥ_{e b} (t) Û_{e} (t - τ) Ĥ_{e b} (τ) Û_{e} (τ - τ^{'}) P Ĥ_{e b} (τ^{'}) Û_{e} (τ^{'} - τ^{''}) Ĥ_{e b} (τ^{''}) \\ - P Ĥ_{e b} (t) Û_{e} (t - τ^{''}) Ĥ_{e b} (τ^{''}) Û_{e} (τ^{''} - τ^{'}) P Ĥ_{e b} (τ^{'}) Û_{e} (τ^{'} - τ) Ĥ_{e b} (τ) \})]) Û_{e}^{†} (t - τ) P {\hat{ρ}}_{e g} (t) . \end{matrix}

(B15)

Taking trace of the above expression over the bath degrees of freedom leads to the following time-local equation:

\frac{d}{d t} {\hat{σ}}_{e g} (t) \approx - \frac{i}{ℏ} Ĥ_{e} {\hat{σ}}_{e g} (t) - (R^{(2)} (t) + R^{(3)} (t) + R^{(4)} (t)) {\hat{σ}}_{e g} (t),

(B16)

where

R^{(2)} (t) = \frac{1}{ℏ^{2}} \int_{0}^{t} d τ {\hat{C}}_{2} (t - τ) Û_{e}^{†} (t - τ),

(B17)

R^{(3)} (t) = - \frac{i}{ℏ^{3}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ {\hat{C}}_{3} (t - τ^{'}, τ^{'} - τ) Û_{e}^{†} (t - τ),

(B18)

\begin{matrix} R^{(4)} (t) & = & - \frac{1}{ℏ^{4}} \int_{0}^{t} d τ^{'} \int_{0}^{τ^{'}} d τ^{''} \int_{0}^{τ^{''}} d τ \\ \times \{{\hat{C}}_{4} (t - τ^{'}, τ^{'} - τ^{''}, τ^{''} - τ)) \\ - {\hat{C}}_{2} (t - τ^{''}) Û_{e} (τ^{''} - τ) {\hat{C}}_{2} (τ - τ^{'}) \\ (- {\hat{C}}_{2} (t - τ^{'}) Û_{e} (τ^{'} - τ) {\hat{C}}_{2} (τ - τ^{''})\} Û_{e}^{†} (t - τ) . \end{matrix}

(B19)

The solution of the above equation provides an alternative method to calculate the line shape, which however requires calculation of ${\hat{σ}}_{e t} (t)$ by direct time evolution. Although this time-local equation approach is not considered in this work, its numerical implementation is feasible for general systems and may prove to be an important alternative way to improve the second order approximation.

APPENDIX C: CALCULATION OF THE 4TH ORDER BATH CORRELATION FUNCTIONS FOR LINEARLY COUPLED HARMONIC OSCILLATOR BATH

For the case of harmonic oscillator bath, the 4th order bath correlation function can be expressed as

\begin{matrix} {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'} + t^{''}) {\hat{B}}_{p_{1} p_{2}} (t^{'} + t^{''}) {\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} \\ = \sum_{n_{1}, n_{2}, n_{3}, n_{4}} ℏ^{4} ω_{n_{1}} ω_{n_{2}} ω_{n_{3}} ω_{n_{4}} \sum_{j_{1}, j_{2}, j_{3}, j_{4}} g_{n_{1} j_{1}} g_{n_{2} j_{2}} g_{n_{3} j_{3}} g_{n_{4} j_{4}} U_{j_{1} p}^{*} U_{j_{1} p_{1}} U_{j_{2} p 1}^{*} U_{j_{2} p_{2}} U_{j_{3} p_{2}}^{*} U_{j_{3} p_{3}} U_{j_{4} p_{3}}^{*} U_{j_{4} q} \\ \times [) {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} e^{- i ω_{n_{1}} (t + t^{'} + t^{''}) - i ω_{n_{2}} (t^{'} + t^{''}) + i ω_{n_{3}} t^{''}} + {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} e^{- i ω_{n_{1}} (t + t^{'} + t^{''}) + i ω_{n_{2}} (t^{'} + t^{''}) - i ω_{n_{3}} t^{''}} \\ + {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} e^{- i ω_{n_{1}} (t + t^{'} + t^{''}) + i ω_{n_{2}} (t^{'} + t^{''}) + i ω_{n_{3}} t^{''}} + {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} e^{i ω_{n_{1}} (t + t^{'} + t^{''}) + i ω_{n_{2}} (t^{'} + t^{''}) - i ω_{n_{3}} t^{''}} \\ + {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} e^{i ω_{n_{1}} (t + t^{'} + t^{''}) - i ω_{n_{2}} (t^{'} + t^{''}) + i ω_{n_{3}} t^{''}} + {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} e^{i ω_{n_{1}} (t + t^{'} + t^{''}) - i ω_{n_{2}} (t^{'} + t^{''}) - i ω_{n_{3}} t^{''}}]) . \end{matrix}

(C1)

Here, each of the trace over the bath can be evaluated as follows:

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ (N_{n_{1}} + 1) (N_{n_{1}} + 2) ⟩ \\ + (1 - δ_{n_{1} n_{2}}) δ_{n_{1} n_{3}} δ_{n_{2} n_{4}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{2}} + 1 ⟩ \\ + (1 - δ_{n_{1} n_{2}}) δ_{n_{1} n_{4}} δ_{n_{2} n_{3}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{2}} + 1 ⟩, \end{matrix}

(C2)

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ {(N_{n_{1}} + 1)}^{2} ⟩ \\ + (1 - δ_{n_{1} n_{3}}) δ_{n_{1} n_{2}} δ_{n_{3} n_{4}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{3}} + 1 ⟩ \\ + (1 - δ_{n_{1} n_{3}}) δ_{n_{1} n_{4}} δ_{n_{2} n_{3}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{3}} ⟩, \end{matrix}

(C3)

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ N_{n_{1}} (N_{n_{1}} + 1) ⟩ \\ + (1 - δ_{n_{1} n_{4}}) δ_{n_{1} n_{2}} δ_{n_{3} n_{4}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{4}} ⟩ \\ + (1 - δ_{n_{1} n_{4}}) δ_{n_{1} n_{3}} δ_{n_{2} n_{4}} ⟨ N_{n_{1}} + 1 ⟩ ⟨ N_{n_{4}} ⟩, \end{matrix}

(C4)

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}}^{†} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ N_{n_{1}} (N_{n_{1}} - 1) ⟩ \\ + (1 - δ_{n_{1} n_{2}}) δ_{n_{1} n_{3}} δ_{n_{2} n_{4}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{2}} ⟩ \\ + (1 - δ_{n_{1} n_{2}}) δ_{n_{1} n_{4}} δ_{n_{2} n_{3}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{2}} ⟩, \end{matrix}

(C5)

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}}^{†} {\hat{b}}_{n_{4}} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ N_{n_{1}}^{2} ⟩ \\ + (1 - δ_{n_{1} n_{3}}) δ_{n_{1} n_{2}} δ_{n_{3} n_{4}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{3}} ⟩ \\ + (1 - δ_{n_{1} n_{3}}) δ_{n_{1} n_{4}} δ_{n_{2} n_{3}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{3}} + 1 ⟩, \end{matrix}

(C6)

\begin{matrix} {T r}_{b} \{{\hat{b}}_{n_{1}}^{†} {\hat{b}}_{n_{2}} {\hat{b}}_{n_{3}} {\hat{b}}_{n_{4}}^{†} {\hat{ρ}}_{b}\} & = & δ_{n_{1} n_{2}} δ_{n_{1} n_{3}} δ_{n_{1} n_{4}} ⟨ N_{n_{1}} (N_{n_{1}} + 1) ⟩ \\ + (1 - δ_{n_{1} n_{4}}) δ_{n_{1} n_{2}} δ_{n_{3} n_{4}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{4}} + 1 ⟩ \\ + (1 - δ_{n_{1} n_{4}}) δ_{n_{1} n_{3}} δ_{n_{2} n_{4}} ⟨ N_{n_{1}} ⟩ ⟨ N_{n_{4}} + 1 ⟩ . \end{matrix}

(C7)

Inserting Eqs. (C2)–(C7) into Eq. (C1) and relabeling the indices of vibrational modes, it is straightforward to derive the following expression:

\begin{matrix} {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'} + t^{''}) {\hat{B}}_{p_{1} p_{2}} (t^{'} + t^{''}) {\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} \\ = \sum_{j_{1}, j_{2}, j_{3}, j_{4}} U_{j_{1} p}^{*} U_{j_{1} p_{1}} U_{j_{2} p_{1}}^{*} U_{j_{2} p_{2}} U_{j_{3} p_{2}}^{*} U_{j_{3} p_{3}} U_{j_{4} p_{3}}^{*} U_{j_{4} q} [) \sum_{n} ℏ^{4} ω_{n}^{4} g_{n j_{1}} g_{n j_{2}} g_{n j_{3}} g_{n j_{4}} \{⟨ N_{n}^{2} + 3 N_{n} + 2 ⟩ e^{- i ω_{n} (t + 2 t^{'} + t^{''})}) \\ + ⟨ N_{n}^{2} + 2 N_{n} + 1 ⟩ e^{- i ω_{n} (t + t^{''})} + ⟨ N_{n}^{2} + N_{n} ⟩ e^{- i ω_{n} (t - t^{''})} + ⟨ N_{n}^{2} - N_{n} ⟩ e^{i ω_{n} (t + 2 t^{'} + t^{''})} + ⟨ N_{n}^{2} ⟩ e^{i ω_{n} (t + t^{''})} + (⟨ N_{n}^{2} + N_{n} ⟩ e^{i ω_{n} (t - t^{''})}\} \\ + \sum_{n \neq n^{'}} ℏ^{4} ω_{n}^{2} ω_{n^{'}}^{2} \{) g_{n j_{1}} g_{n^{'} j_{2}} g_{n j_{3}} g_{n^{'} j_{4}} () ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{- i ω_{n} (t + t^{'}) - i ω_{n^{'}} (t^{'} + t^{''})} + ⟨ N_{n} ⟩ ⟨ N_{n^{'}} ⟩ e^{i ω_{n} (t + t^{'}) + i ω_{n^{'}} (t^{'} + t^{''})} + ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} ⟩ e^{- i ω_{n} (t + t^{'}) + i ω_{n^{'}} (t^{'} + t^{''})} \\ + ⟨ N_{n} ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{i ω_{n} (t + t^{'}) - i ω_{n^{'}} (t^{'} + t^{''})})) + g_{n j_{1}} g_{n^{'} j_{2}} g_{n^{'} j_{3}} g_{n j_{4}} () ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{- i ω_{n} (t + t^{'} + t^{''}) - i ω_{n^{'}} t^{'}} + ⟨ N_{n} ⟩ ⟨ N_{n}^{'} ⟩ e^{i ω_{n} (t + t^{'} + t^{''}) + i ω_{n^{'}} t^{'}} \\ + ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} ⟩ e^{- i ω_{n} (t + t^{'} + t^{''}) + i ω_{n^{'}} t^{'}} + ⟨ N_{n} ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{i ω_{n} (t + t^{'} + t^{''}) - i ω_{n^{'}} t^{'}})) + g_{n j_{1}} g_{n j_{2}} g_{n^{'} j_{3}} g_{n^{'} j_{4}} () ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{- i ω_{n} t - i ω_{n^{'}} t^{''}} \\ + ⟨ N_{n} ⟩ ⟨ N_{n^{'}} ⟩ e^{i ω_{n} t + i ω_{n^{'}} t^{''}} + ⟨ N_{n} + 1 ⟩ ⟨ N_{n^{'}} ⟩ e^{- i ω_{n} t + i ω_{n^{'}} t^{''}} (+ ⟨ N_{n} ⟩ ⟨ N_{n^{'}} + 1 ⟩ e^{i ω_{n} t - i ω_{n^{'}} t^{''}}))\}]) . \end{matrix}

(C8)

In the above expression, $⟨ N_{n} ⟩ = 1 / (e^{β ℏ ω_{n}} - 1)$ and $⟨ N_{n}^{2} ⟩ = 1 / (e^{β ℏ ω_{n}} - 1) + 2 / {(e^{β ℏ ω_{n}} - 1)}^{2} = 2 {⟨ N_{n} ⟩}^{2} + ⟨ N_{n} ⟩$ ⁠. Using these expressions, it is straightforward to show that the first sum over n in Eq. (C8) corresponds to the case of the second sum if n′ were the same as n. Thus, the first sum over n can be incorporated into the second sum by removing the restriction, n ≠ n′. Employing this expression and the one used for calculating Eq. (52), the fourth order term in Eq. (42) can be expressed as

\begin{matrix} {T r}_{b} \{{\hat{B}}_{p p_{1}} (t + t^{'} + t^{''}) {\hat{B}}_{p_{1} p_{2}} (t^{'} + t^{''}) {\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} - {T r}_{b} \{{\hat{B}}_{p p_{1}} (t) {\hat{B}}_{p_{1} p_{2}} {\hat{ρ}}_{b}\} {T r}_{b} \{{\hat{B}}_{p_{2} p_{3}} (t^{''}) {\hat{B}}_{p_{3} q} {\hat{ρ}}_{b}\} \\ = \frac{ℏ^{2}}{π^{2}} \sum_{j_{1}, j_{2}, j_{3}, j_{4}} U_{j_{1} p}^{*} U_{j_{1} p_{1}} U_{j_{2} p_{1}}^{*} U_{j_{2} p_{2}} U_{j_{3} p_{2}}^{*} U_{j_{3} p_{3}} U_{j_{4} p_{3}}^{*} U_{j_{4} q} \\ \times [\int_{0}^{\infty} d ω J_{j_{1} j_{3}} (ω) \int_{0}^{\infty} d ω^{'} J_{j_{2} j_{4}} (ω^{'})) \{\frac{1}{(1 - e^{- β ℏ ω})} \frac{1}{(1 - e^{- β ℏ ω^{'}})} e^{- i ω t - i (ω + ω^{'}) t^{'} - i ω^{'} t^{''}}) \\ + \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{e^{- β ℏ ω^{'}}}{(1 - e^{- β ℏ ω^{'}})} e^{i ω t + i (ω + ω^{'}) t^{'} + i ω^{'} t^{''}} + \frac{1}{(1 - e^{- β ℏ ω})} \frac{e^{- β ℏ ω^{'}}}{(1 - e^{- β ℏ ω^{'}})} e^{- i ω t - i (ω - ω^{'}) t^{'} + i ω^{'} t^{''}} \\ (+ \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{1}{(1 - e^{- β ℏ ω^{'}})} e^{i ω t + i (ω - ω^{'}) t^{'} - i ω^{'} t^{''}}\} + \int_{0}^{\infty} d ω J_{j_{1} j_{4}} (ω) \int_{0}^{\infty} d ω^{'} J_{j_{2} j_{3}} (ω^{'}) \\ \times \{\frac{1}{(1 - e^{- β ℏ ω})} \frac{1}{(1 - e^{- β ℏ ω^{'}})} e^{- i ω t - i (ω + ω^{'}) t^{'} - i ω t^{''}}) + \frac{e^{- β ℏ ω}}{(1 - e^{- β ℏ ω})} \frac{e^{- β ℏ ω^{'}}}{(1 - e^{- β ℏ ω^{'}})} e^{i ω t + i (ω + ω^{'}) t^{'} + i ω t^{''}} \\ (+ \frac{1}{(1 - e^{- β ℏ ω})} \frac{e^{- β ℏ ω^{'}}}{(1 - e^{- β ℏ ω^{'}})} e^{- i ω t - i (ω - ω^{'}) t^{'} - i ω t^{''}} (+ \frac{e^{- β ℏ ω}}{(1 - e^{β ℏ ω})} \frac{1}{(1 - e^{- β ℏ ω^{'}})} e^{i ω t + i (ω - ω^{'}) t^{'} + i ω t^{''}}\}] . \end{matrix}

(C9)

REFERENCES

1.

J. R. G.

Breene

,

Theories of Spectral Line Shape

(

John Wiley & Sons, Inc.

,

New York

,

1981

).

Google Scholar

Crossref

2.

S.

Mukamel

,

Principles of Nonlinear Spectroscopy

(

Oxford University Press

,

New York

,

1995

).

Google Scholar

3.

K.

Blum

,

Density Matrix Theory and Applications

(

Plenum Press

,

New York

,

1996

).

Google Scholar

Crossref

4.

I.

de Vega

and

D.

Alonso

,

Rev. Mod. Phys.

89

,

015001

(

2017

).

https://doi.org/10.1103/revmodphys.89.015001

Google Scholar

Crossref

5.

S. J.

Jang

and

B.

Mennucci

,

Rev. Mod. Phys.

90

,

035003

(

2018

).

https://doi.org/10.1103/revmodphys.90.035003

Google Scholar

Crossref

6.

A. S.

Davydov

,

Theory of Molecular Excitons

(

Plenum Press

,

New York, London

,

1971

).

Google Scholar

Crossref

7.

I. I.

Abram

and

R.

Silbey

,

J. Chem. Phys.

63

,

2317

(

1975

).

https://doi.org/10.1063/1.431683

Google Scholar

Crossref

8.

M.

Blume

,

Phys. Rev.

174

,

351

(

1968

).

https://doi.org/10.1103/physrev.174.351

Google Scholar

Crossref

9.

V. M.

Agranovich

and

E. R. M.

Hochstrasser

,

Spectroscopy and Excitation Dynamics of Condensed Molecular Systems

, Modern Problems in Condensed Matter Sciences Vol. 4 (

North-Holland

,

Amsterdam

,

1982

).

Google Scholar

10.

Organic Molecular Aggregates

, edited by

P.

Reineker

,

H.

Haken

, and

H. C.

Wolf

(

Springer-Verlag

,

Berlin

,

1983

).

Google Scholar

Crossref

11.

H.

Fidder

,

J.

Knoester

, and

D. A.

Wiersma

,

J. Chem. Phys.

95

,

7880

(

1991

).

https://doi.org/10.1063/1.461317

Google Scholar

Crossref

12.

N. J.

Hestand

and

F. C.

Spano

,

Chem. Rev.

118

,

7069

(

2018

).

https://doi.org/10.1021/acs.chemrev.7b00581

Google Scholar

Crossref

PubMed

13.

M.

Hoffmann

and

Z. G.

Soos

,

Phys. Rev. B

66

,

024305

(

2002

).

https://doi.org/10.1103/physrevb.66.024305

Google Scholar

Crossref

14.

F. C.

Spano

,

Annu. Rev. Phys. Chem.

57

,

217

(

2006

).

https://doi.org/10.1146/annurev.physchem.57.032905.104557

Google Scholar

Crossref

PubMed

15.

M.

Cho

,

Chem. Rev.

108

,

1331

(

2008

).

https://doi.org/10.1021/cr078377b

Google Scholar

Crossref

PubMed

16.

V.

Butkus

,

D.

Abramavicius

,

A.

Gelzinis

, and

L.

Valkunas

,

Lith. J. Phys.

50

,

267

(

2010

).

https://doi.org/10.3952/lithjphys.50305

Google Scholar

Crossref

17.

A. M.

van Oijen

et al,

Science

285

,

400

(

1999

).

https://doi.org/10.1126/science.285.5426.400

Google Scholar

Crossref

PubMed

18.

R.

Hildner

,

U.

Lemmer

,

U.

Scherf

,

M.

van Heel

, and

J.

Köhler

,

Adv. Mater.

19

,

1978

(

2007

).

https://doi.org/10.1002/adma.200602718

Google Scholar

Crossref

19.

C.

Hofmann

et al,

Phys. Rev. Lett.

90

,

013004

(

2003

).

https://doi.org/10.1103/physrevlett.90.013004

Google Scholar

Crossref

PubMed

20.

J.

Hofkens

et al,

J. Am. Chem. Soc.

122

,

9278

(

2000

).

https://doi.org/10.1021/ja0012570

Google Scholar

Crossref

21.

S.

Jang

,

M. D.

Newton

, and

R. J.

Silbey

,

Phys. Rev. Lett.

92

,

218301

(

2004

).

https://doi.org/10.1103/physrevlett.92.218301

Google Scholar

Crossref

PubMed

22.

S.

Jang

,

S.

Hoyer

,

G. R.

Fleming

, and

K. B.

Whaley

,

Phys. Rev. Lett.

113

,

188102

(

2014

).

https://doi.org/10.1103/physrevlett.113.188102

Google Scholar

Crossref

PubMed

23.

K.

Faid

and

R. F.

Fox

,

J. Chem. Phys.

88

,

4579

(

1988

).

https://doi.org/10.1063/1.453770

Google Scholar

Crossref

24.

M.

Yang

,

J. Chem. Phys.

123

,

124705

(

2005

).

https://doi.org/10.1063/1.2046668

Google Scholar

Crossref

PubMed

25.

A.

Gelzinis

,

D.

Abramavicius

, and

L.

Valkunas

,

J. Chem. Phys.

142

,

154107

(

2015

).

https://doi.org/10.1063/1.4918343

Google Scholar

Crossref

PubMed

26.

L.

Chen

,

R.

Zheng

,

Q.

Shi

, and

Y.

Yan

,

J. Chem. Phys.

131

,

094502

(

2009

).

https://doi.org/10.1063/1.3213013

Google Scholar

Crossref

PubMed

27.

M.

Schröder

,

U.

Kleinekathöfer

, and

M.

Schreiber

,

J. Chem. Phys.

124

,

084903

(

2006

).

https://doi.org/10.1063/1.2171188

Google Scholar

Crossref

PubMed

28.

T.-C.

Dinh

and

T.

Renger

,

J. Chem. Phys.

142

,

034104

(

2015

).

https://doi.org/10.1063/1.4904928

Google Scholar

Crossref

PubMed

29.

T.-C.

Dinh

and

T.

Renger

,

J. Chem. Phys.

145

,

034105

(

2016

).

https://doi.org/10.1063/1.4958322

Google Scholar

Crossref

PubMed

30.

P.

Kumar

and

S.

Jang

,

J. Chem. Phys.

138

,

135101

(

2013

).

https://doi.org/10.1063/1.4795824

Google Scholar

Crossref

PubMed

31.

U.

Kleinekathöfer

,

J. Chem. Phys.

121

,

2505

(

2004

).

https://doi.org/10.1063/1.1770619

Google Scholar

Crossref

PubMed

32.

A.

Pomyalov

and

D. J.

Tanor

,

J. Chem. Phys.

123

,

204111

(

2005

).

https://doi.org/10.1063/1.2121649

Google Scholar

Crossref

PubMed

33.

A.

Montoya-Castillo

,

T. C.

Berkelbach

, and

D. R.

Reichman

,

J. Chem. Phys.

143

,

194108

(

2015

).

https://doi.org/10.1063/1.4935443

Google Scholar

Crossref

PubMed

34.

J. H.

Fetherolf

and

T. C.

Berkelbach

,

J. Chem. Phys.

147

,

244109

(

2017

).

https://doi.org/10.1063/1.5006824

Google Scholar

Crossref

PubMed

35.

S.

Jang

,

J.

Cao

, and

R. J.

Silbey

,

J. Chem. Phys.

116

,

2705

(

2002

).

https://doi.org/10.1063/1.1445105

Google Scholar

Crossref

36.

B. B.

Laird

,

J.

Budimir

, and

J. L.

Skinner

,

J. Chem. Phys.

94

,

4391

(

1991

).

https://doi.org/10.1063/1.460626

Google Scholar

Crossref

37.

B. B.

Laird

and

J. L.

Skinner

,

J. Chem. Phys.

94

,

4405

(

1991

).

https://doi.org/10.1063/1.460627

Google Scholar

Crossref

38.

H. P.

Breuer

,

B.

Kappler

, and

F.

Petruccione

,

Ann. Phys.

291

,

36

(

2001

).

https://doi.org/10.1006/aphy.2001.6152

Google Scholar

Crossref

39.

A. A.

Golosov

and

D. R.

Reichman

,

Chem. Phys.

296

,

129

(

2004

).

https://doi.org/10.1016/j.chemphys.2003.10.009

Google Scholar

Crossref

40.

Q.

Shi

and

E.

Geva

,

J. Chem. Phys.

119

,

12063

(

2003

).

https://doi.org/10.1063/1.1624830

Google Scholar

Crossref

41.

M.-L.

Zhang

,

B. J.

Ka

, and

E.

Geva

,

J. Chem. Phys.

125

,

044106

(

2006

).

https://doi.org/10.1063/1.2218342

Google Scholar

Crossref

42.

R. X.

Xu

and

Y. J.

Yan

,

Phys. Rev. E

75

,

031107

(

2007

).

https://doi.org/10.1103/physreve.75.031107

Google Scholar

Crossref

43.

M.

Schroder

,

M.

Schreiber

, and

U.

Kleinekathofer

,

J. Chem. Phys.

126

,

114102

(

2007

).

https://doi.org/10.1063/1.2538754

Google Scholar

Crossref

PubMed

44.

R.

Doll

et al,

Chem. Phys.

347

,

243

(

2008

).

https://doi.org/10.1016/j.chemphys.2007.09.003

Google Scholar

Crossref

45.

H. T.

Chen

,

T. C.

Berkelbach

, and

D. R.

Reichman

,

J. Chem. Phys.

144

,

154106

(

2016

).

https://doi.org/10.1063/1.4946809

Google Scholar

Crossref

PubMed

46.

S.

Jang

and

R. J.

Silbey

,

J. Chem. Phys.

118

,

9312

(

2003

).

https://doi.org/10.1063/1.1569239

Google Scholar

Crossref

47.

A. A.

Golosov

and

D. R.

Reichman

,

J. Chem. Phys.

115

,

9848

(

2001

).

https://doi.org/10.1063/1.1412611

Google Scholar

Crossref

48.

J. L.

Wu

and

J.

Cao

,

J. Chem. Phys.

139

,

044102

(

2013

).

https://doi.org/10.1063/1.4812781

Google Scholar

Crossref

PubMed

49.

M. G.

Mavros

and

T.

van Voorhis

,

J. Chem. Phys.

141

,

054112

(

2014

).

https://doi.org/10.1063/1.4891669

Google Scholar

Crossref

PubMed

50.

M.

Xu

,

Y.

Yan

,

Y.

Liu

, and

Q.

Shi

,

J. Chem. Phys.

148

,

164101

(

2018

).

https://doi.org/10.1063/1.5022761

Google Scholar

Crossref

PubMed

51.

A. A.

Golosov

and

D. R.

Reichman

,

J. Chem. Phys.

115

,

9862

(

2001

).

https://doi.org/10.1063/1.1412612

Google Scholar

Crossref

52.

G.

Ritschel

and

A.

Eisfeld

,

J. Chem. Phys.

141

,

094101

(

2014

).

https://doi.org/10.1063/1.4893931

Google Scholar

Crossref

PubMed

2019

Author(s)

Fourth order expressions for the electronic absorption lineshape of molecular excitons

I. INTRODUCTION

II. GENERAL FORMALISM

III. LINE SHAPE EXPRESSIONS FOR A LINEARLY COUPLED HARMONIC OSCILLATOR BATH MODEL

IV. MODEL CALCULATIONS

V. CONCLUSION

ACKNOWLEDGMENTS

APPENDIX A: SECOND ORDER APPROXIMATION FOR THE EXCITON POPULATION

APPENDIX B: EXACT FORMAL EXPRESSION AND FOURTH ORDER APPROXIMATION FOR THE TIME EVOLUTION OF PROJECTED DENSITY OPERATOR

APPENDIX C: CALCULATION OF THE 4TH ORDER BATH CORRELATION FUNCTIONS FOR LINEARLY COUPLED HARMONIC OSCILLATOR BATH

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Fourth order expressions for the electronic absorption lineshape of molecular excitons

I. INTRODUCTION

II. GENERAL FORMALISM

III. LINE SHAPE EXPRESSIONS FOR A LINEARLY COUPLED HARMONIC OSCILLATOR BATH MODEL

IV. MODEL CALCULATIONS

V. CONCLUSION

ACKNOWLEDGMENTS

APPENDIX A: SECOND ORDER APPROXIMATION FOR THE EXCITON POPULATION

APPENDIX B: EXACT FORMAL EXPRESSION AND FOURTH ORDER APPROXIMATION FOR THE TIME EVOLUTION OF PROJECTED DENSITY OPERATOR

APPENDIX C: CALCULATION OF THE 4TH ORDER BATH CORRELATION FUNCTIONS FOR LINEARLY COUPLED HARMONIC OSCILLATOR BATH

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only