Reduced quantum dynamics with arbitrary bath spectral densities: hierarchical equations of motion based on several different bath decomposition schemes.

Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes Hao Liu, Lili Zhu, Shuming Bai, and Qiang Shi Citation: The Journal of Chemical Physics 140, 134106 (2014); doi: 10.1063/1.4870035 View online: http://dx.doi.org/10.1063/1.4870035 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/13?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Non-Markovian reduced dynamics based upon a hierarchical effective-mode representation J. Chem. Phys. 137, 144107 (2012); 10.1063/1.4752078 Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems J. Chem. Phys. 134, 244106 (2011); 10.1063/1.3602466 Effective-mode representation of non-Markovian dynamics: A hierarchical approximation of the spectral density. II. Application to environment-induced nonadiabatic dynamics J. Chem. Phys. 131, 124108 (2009); 10.1063/1.3226343 Spectral and entropic characterizations of Wigner functions: Applications to model vibrational systems J. Chem. Phys. 129, 094103 (2008); 10.1063/1.2968607 Dynamics of quantum dissipation systems interacting with fermion and boson grand canonical bath ensembles: Hierarchical equations of motion approach J. Chem. Phys. 126, 134113 (2007); 10.1063/1.2713104

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.12.65.247 On: Sun, 27 Apr 2014 18:03:43

THE JOURNAL OF CHEMICAL PHYSICS 140, 134106 (2014)

Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes Hao Liu, Lili Zhu, Shuming Bai, and Qiang Shia) Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100190, China

(Received 5 November 2013; accepted 18 March 2014; published online 4 April 2014) We investigated applications of the hierarchical equation of motion (HEOM) method to perform high order perturbation calculations of reduced quantum dynamics for a harmonic bath with arbitrary spectral densities. Three different schemes are used to decompose the bath spectral density into analytical forms that are suitable to the HEOM treatment: (1) The multiple Lorentzian mode model that can be obtained by numerically fitting the model spectral density. (2) The combined Debye and oscillatory Debye modes model that can be constructed by fitting the corresponding classical bath correlation function. (3) A new method that uses undamped harmonic oscillator modes explicitly in the HEOM formalism. Methods to extract system-bath correlations were investigated for the above bath decomposition schemes. We also show that HEOM in the undamped harmonic oscillator modes can give detailed information on the partial Wigner transform of the total density operator. Theoretical analysis and numerical simulations of the spin-Boson dynamics and the absorption line shape of molecular dimers show that the HEOM formalism for high order perturbations can serve as an important tool in studying the quantum dissipative dynamics in the intermediate coupling regime. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4870035] I. INTRODUCTION

A large number of problems in chemical dynamics involve studying the dynamics of a quantum system coupled to a complex bath.1 Typical examples include charge and energy transfer processes,2–4 linear and nonlinear spectroscopy,5, 6 and charge transport.7 In literature, several powerful approximate methods based on the Fermi’s golden rule (FGR), for example, the Marcus theory8, 9 and Förster theory,10 have been applied in such problems. At the same time, the reduced quantum dynamics methods such as the Redfield theory,11, 12 and various different forms of the second order generalized quantum master equations5, 13 (GQMEs) are also widely employed. Both of these methods are based on second order perturbation: FGR treats the inter-system couplings to the second order, while the Redfield theory and the second order GQMEs apply the second order perturbation for the system-bath coupling, with or without further Markovian approximations. A challenging problem for the above theories is the quantum dissipative dynamics in the so-called intermediate coupling regime, where the inter-system couplings and the system-bath couplings are of similar strength.3, 4 It has been shown that the above mentioned second order perturbation methods become invalid in calculations of excitation energy transfer rates,14, 15 as well as linear and nonlinear spectroscopy.16–18 Different schemes to go beyond the second order approximations are proposed in literature. A natural way is to dea) [email protected]

0021-9606/2014/140(13)/134106/11/$30.00

rive reduced equations of motion in higher order perturbations. For example, the fourth order GQMEs have been derived by several groups.19–21 However, the resulted fourth order equations of motion become considerably complicated, and their extensions to higher orders are difficult. Another important problem is that, there is no systematic test of whether the fourth order perturbation is adequate when the second order approximation becomes invalid. In fact, when the second order perturbation fails, convergence of the fourth order approximations should also be checked to find out whether even higher order perturbations are needed.17, 22–24 A special form of non-perturbative quantum master equation has also been proposed based on the Padé approximation, where the second and fourth order expansions are used to construct an approximate infinite order perturbation.25, 26 These methods were found to work well when the collective bath dynamics is close to Gaussian-Markovian, while their applicability in general cases depends on the accuracy of the Padé approximation with only the second and fourth order terms. An alternative approach is to explicitly include a collective bath coordinate into the system degrees of freedom, such that the weak coupling approximation could be applied in the newly defined system-bath problem with the extended quantum subsystem.12, 27 More recently, the polaron transformation approach has also been used to reduce the effective system-bath coupling, such that a new set of perturbative quantum master equations with a larger range of applicability can be derived.28–30 In recent years, the Liouville space hierarchical equation of motion (HEOM) method31–36 has emerged as a

140, 134106-1

© 2014 AIP Publishing LLC


134106-2

Liu et al.

powerful tool in simulating dissipative quantum dynamics in condensed phase. It has been realized that the truncation at the Lth order of hierarchical is equivalent to N = 2L order perturbation theory.17, 37 The HEOM formalism can thus provide a convenient way to perform high-order perturbations in the system-bath coupling. A major advantage of applying the HEOM formalism for problems in the intermediate coupling regime is that, since the system-bath coupling is comparable but not much larger than the inter-system coupling, the equations of motion may not need very high orders to converge. At the same time, the HEOM formalism provides a systematic way to check the convergence with perturbation orders, which is not available in previous high order perturbation approximations. Up to now, the HEOM approach has been mostly applied to problems with the Debye and Brownian oscillator (Lorentzian) spectral densities, or simple combinations of a few of them, due to the reason that the bath correlation function can be written analytically into a sum of exponential functions. In the case of general spectral densities, Meier and Tannor38 have proposed a scheme to use a sum of Lorentzian modes to numerically fit a given spectral density, and this approach has been applied successfully to second order GQMEs.27, 38–40 The same decomposition scheme can certainly be used to treat general spectral densities in the HEOM method. Later, Kleinekathöfer and co-workers have used combined Debye and oscillatory Debye modes to fit the numerically calculated classical bath correlation functions, which results in a similar but different decomposition scheme.41 More recently, Kramer and co-workers have applied the high temperature approximation of the HEOM for simulations of energy transfer dynamics with such spectral densities.42 In this work, we will use the HEOM combined with the above decomposition schemes for arbitrary bath spectral densities to investigate the convergence of high order perturbations. Besides using the above mentioned dissipative models (multiple Lorentzian modes, combinations of Debye and oscillatory Debye modes) to decompose the spectral density, it is often desirable to include some discrete harmonic oscillator modes in the system-bath coupling to describe electronicvibrational coupling to specific vibrational modes. To treat these discrete modes, a simple way is to include them in the system degrees of freedom, resulting in an extended system Hamiltonian. In this study, we show that such undamped harmonic oscillator modes can also be treated effectively within the HEOM formalism, and the Wigner distribution47 of these modes can be extracted from the auxiliary density operators (ADOs). The remaining sections of this paper are arranged as follows. In Sec. II, the model Hamiltonian, HEOM for an arbitrary spectral density based on the multiple Lorentzian mode decomposition scheme,38 and the combined Debye and oscillatory Debye modes41 are first presented. The HEOM formalism with undamped harmonic oscillator modes, and the relations between the corresponding ADOs with the partial Wigner transform of the harmonic oscillator modes are also derived. In Sec. III, numerical simulations are presented to calculate the spin-Boson dynamics and the absorption spectra

J. Chem. Phys. 140, 134106 (2014)

of model molecular dimers with different bath spectral densities. Conclusions and discussions are made in Sec. IV. II. THEORY

We consider a quantum system coupled to a harmonic bath.43–45 The total system-bath Hamiltonian can be written as ⎡ 2 ⎤ 2 pˆ j2 ˆ c p 1 j ⎣ ⎦, ˆ + ˆ Hˆ = + Vs (q) + ωj2 xˆj − 2 f (q) 2 2 2 ω j j (1) where qˆ and pˆ describe the system degree of freedom; xˆj and pˆ j are the mass-weighted coordinate and momentum of the jth bath mode with frequency ωj , cj is the coupling coefficient ˆ and the jth bath mode cobetween the system operator f (q) ordinate xˆj . The Hamiltonian in Eq. (1) can be rewritten as Hˆ = Hˆ s + Hˆ b + Hˆ sb + Hˆ ren ,

(2)

pˆ 2 ˆ + Vs (q) Hˆ s = 2

(3)

where

is the system Hamiltonian, Hˆ b =

pˆ j2

1 + ωj2 xˆj2 2 2

j

is the bath Hamiltonian, ˆ Hˆ sb = − cj xˆj ⊗ qˆ = −Fˆ ⊗ f (q)

(4)

(5)

j

describes the system-bath interaction, and Hˆ ren =

cj2 j

2ωj2

ˆ 2 f (q)

(6)

is the renormalization term. The system-bath interaction can be characterized by the spectral density J(ω) defined as J (ω) =

2 π cj δ(ω − ωj ). 2 j ωj

(7)

The time correlation function C(t) for the collective bath co ordinate F = j cj xˆj is related to the spectral density via the fluctuation-dissipation theorem, and can be expanded into an exponential-like series:

¯ ∞ e−iωt C(t) = Fˆ (t)Fˆ B = dωJ (ω) π −∞ 1 − e−βω = dk e−νk t , for t > 0. (8) k

While for t < 0, a different expansion scheme should be used, and it can be proven that C(−t) is the complex conjugate of C(t).


134106-3

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

Equation (8) is an important step in deriving the HEOM formalism. When using dissipative modes with the Debye or Lorentzian spectral density, it can be shown that the expansion is exact when taking into infinite number of expansion terms. In real calculations, only a finite number of expansion terms can be treated, terms with large ν k values can be treated using the Markovian approximation,34, 36 which is found to be valid when the absolute value of ν k is much faster than the system time scale.36 The general form of HEOM can then be written as35 d i ˆ nk νk ρn − ρn+k ρn = − iL + f (q), dt ¯ k k i ˆ n−k − dk∗ ρn−k f (q)), ˆ − nk (dk f (q)ρ ¯ k

(9)

where Lρˆ = ¯−1 [Hˆ s + Hˆ ren , ρ], ˆ ρ n s are generalized density operators of the system degree of freedom, where − the subscripts n, n+ k , nk denote {n1 , . . . , nk , . . . }, {n1 , . . . , nk + 1, . . . }, and {n1 , . . . , nk − 1, . . . }, respectively. The density operator ρ 0 = ρ 0, ..0, .., 0 corresponds to the reduced density matrix of the system, others are the auxiliary density operators. As stated above, the N = 2L order perturbation can be obtained by truncation of the HEOM at the Lth order (L = k nk ) of hierarchy.17, 37 The main issue to apply the HEOM formalism for general spectral densities is then how to obtain the exponential expansion of the bath correlation function in Eq. (8), and how to solve the resulted Eq. (9) efficiently. Generally, a first step is to decompose the spectral density into suitable analytically forms, where the exponential expansion in Eq. (8) can be calculated easily. We will discuss two commonly used decomposition schemes in Secs. II A and II B, and also a scheme to incorporate undamped harmonic oscillator modes within the HEOM formalism in Sec. II C. For the reason of completeness, detailed expressions of the HEOM for the three bath decomposition schemes are given in Appendixes A–C. Sections II D and II E deal with the calculation of systembath correlations from the ADOs. In a previous paper,46 we have shown that higher moments of the collective bath coordinate can be calculated using the ADOs, while the distribution function of collective bath coordinate was calculated numerically from these high order moments. Here, we go a further step by showing that such distribution function can be calculated directly from the ADOs, which is presented in Sec. II D. In the case of undamped harmonic oscillator modes, we show in Sec. II E that the partial Wigner transform47 with respect to the harmonic oscillator modes can also be calculated directly using the ADOs.

The idea of combining this decomposition scheme with the HEOM formalism has been discussed previously.15, 36, 48 Details of the HEOM in this bath decomposition scheme are given in Appendix A. B. Sum of Debye and oscillatory Debye modes

Another type of decomposition scheme has recently been employed by Kleinekathöfer and co-workers, which was obtained by fitting the classical correlation functions to a summation of exponential and oscillatory exponential decays.41 In this scheme, the bath spectral density is written as J (ω) =

Na ηa,j γa,j ω 2 ω2 + γa,j j =1

+

Nb

j =1

ηb,j γb,j ω ηb,j γb,j ω + 2 2 2 (ω − ωb,j ) + γb,j (ω + ωb,j )2 + γb,j

the classical bath correlation function of which is given by a summation of exponential and oscillatory exponential decays: Ccl (t) =

Na

ηa,j e−γa,j t +

j =1

J (ω) =

N

ω . (10) pj 2 2 (ω + j ) + j (ω − j )2 + j2 j =1

Nb

ηb,j cos(ωb,j t)e−γb,j t .

(12)

j =1

We note that this form of correlation functions has also been used to fit the autocorrelation functions from experiments.49, 50 The high temperature approximation of HEOM with the decomposition scheme in Eq. (11) has been given in Ref. 42, while the full quantum HEOM is given in Appendix B. C. Using undamped harmonic oscillator modes

In some problems, such as when considering an individual harmonic oscillator mode that couples strongly to the system degree of freedom, it is necessary to include discrete harmonic oscillator modes in the spectral density. Although a natural way is to redefine an extended system Hamiltonian to include these modes, we show that the HEOM formalism is also capable to accommodate such undamped harmonic oscillator modes explicitly. In such case, the exponential constants ν k s in Eq. (8) are purely imaginary, and there are no Matsubara terms. Assuming N undamped harmonic oscillator modes, the bath correlation function in Eq. (8) can be written as C(t) =

N ¯cj2 j =1

Meier and Tannor have proposed a decomposition scheme where the spectral density J(ω) is numerically fitted into a sum of (under damped) Lorentzians:

,

(11)

A. Sum of Lorentzian modes 38

+

4ωj

[coth(β¯ωj /2) − 1]eiωj t

N ¯cj2 j =1

4ωj

[coth(β¯ωj /2) + 1]e−iωj t ,

(13)

which results in the HEOM presented in Appendix C. We will show later in Sec. II E that the HEOM formulated using undamped modes actually gives a full description of the systembath correlations.


134106-4

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

D. Extracting the bath information: General cases

As found in previous studies, the ADOs in the HEOM formalism contain useful information about the system-bath correlations.15, 17, 24, 51 Quantitative ways to calculate such system-bath correlations were also found: (1) In the high temperature limit with the Debye spectral density, the ADOs in the HEOM can be explained as expansion coefficients of the Wigner distribution function of the overdamped collective bath coordinate in a predefined basis set.36 (2) In the case of general spectral densities and arbitrary temperatures, we have shown that higher order expectation values of the collective bath coordinate can be calculated with the ADOs.46 In Ref. 46, we have also constructed the distribution of the collective bath coordinate from the higher order moments using numerical fitting. In this work, we will present a new analytical equation that can be used to calculate the distribution function of the collective bath coordinate directly. Using the definition of the collective bath coordinate Fˆ defined in Eq. (5), the higher order expectation values of the Fˆ are defined as X(n) = TrB [Fˆ n ρT (t)].

(14)

Sec. II C, the partial Wigner distribution of the total density operator can be calculated from the ADOs in Eq. (C1). The partial Wigner distribution of the total density operator is defined as

1 i T de− ¯ p · x + /2|ρˆT |x − /2, ρw (x, p) = (2π ¯)N (19) where N is the total number of degrees of freedom, ρ T is the density operator of the whole system. The equation of motion for ρw (t) can be obtained by doing the partial Wigner transform of the quantum Liouville equation.35, 52 In the case of linear coupling to the Harmonic bath as described by the Hamiltonian in Eq. (1), the following equations of motion can be derived:53–55 ∂ i i ˆ ρw ] ρw (t) = − [Hs , ρw ] + Lj ρw + cj xj [f (q), ∂t ¯ ¯ j j −

cj ∂ {f (q), ˆ ρw } , 2 ∂pj j

(20)

where

(n)

We note that X is an operator of the system degree of freedom. We then show that a “distribution function” of collective coordinate Fˆ can be constructed from the high order moments calculated in Eq. (14). To this end, the distribution is expanded with the following basis set n (X): an n (X), (15) P (X) =

Lj = −pj

n (X) =

1 2π

1/2 √

1

2 −X Hn (X/ 2A0 )e 2A0 ,

2n n! where Hn (x) is the nth order Hermite polynomial, and

1 ∞ β¯ω . A0 = dωJ (ω) coth π 0 2

L= (16)

(17)

The expansion coefficients an can be calculated from the ADOs by n! (−1)n ρ an = √ , nk = n. (18) n nk ! n!A0 n th

(21)

The eigenfunctions of Lj can be obtained as follows (for simplicity, we drop the subscript j):

n

∂ ∂ + ωj2 xj . ∂xj ∂pj

∂ ∂ 1 L = −p˜ + x˜ , ω ∂ x˜ ∂ p˜

(22)

where the reduced coordinates and momenta x˜ and p˜ are defined as 2ω 2 , p˜ = p . x˜ = x ¯ coth(¯ωβ/2) ¯ω coth(¯ωβ/2) (23) By making use of the following transformation, 2 2 p˜ 2 x˜ p˜ 2 x˜ + L exp − − , L˜ = exp 4 4 4 4

(24)

n order

The detailed derivation of Eq. (18) is presented in Appendix D. Equations (15)–(18) extend the findings in a previous paper,46 in that we can now calculate the distribution function P(X) directly, without the need to calculate the higher order moments in Eq. (14) as a first step. Another advantage is that we can employ a different scaling factor in the previously developed on-the-fly filtering algorithm36 by defining 1 nk − 2 ρn . This new renormalization factor is ρ˜n = k nk !A0 found to be advantageous in some applications.

it is found that,

E. Undamped harmonic oscillator modes: Relation to the partial Wigner transform

The transformed operator L˜ is then diagonalized to obtain

In this subsection, we show that, in the equation of motion with the undamped harmonic oscillator modes derived in

L˜ = −ab+ + a + b,

(25)

where a=

x˜ ∂ + , ∂ x˜ 2

a+ = −

x˜ ∂ + , ∂ x˜ 2

(26a)

b=

∂ p˜ + , ∂ p˜ 2

b+ = −

∂ p˜ + . ∂ p˜ 2

(26b)

L˜ = ic1+ c1 − ic2+ c2 ,

(27)


134106-5

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

where the creation and annihilation operators of the two eigenmodes c1± and c2± are given by 1 c1 = √ (−ib + a), 2 1 c2 = √ (b − ia), 2

1 c1+ = √ (ib+ + a + ), 2 c2+

in Eq. (C1) via ρ˜m,n =

mj nj

(−1) i

¯cj2 coth(¯ωj β/2) 4ωj

j

(28)

1 = √ (b+ + ia + ). 2

−(mj +nj )/2 ρm,n ,

(34) Eqs. (31) and (34) thus provide a way to calculate the partial Wigner transform ρw from the ADOs.

The eigenfunctions of L are given by φm,n = √

III. RESULTS

1 m!n!

(c1+ )m (c2+ )n φ00

,

(29)

√ 2 2 where φ00 = 1/ 2π e(−x˜ /4−p˜ /4) denotes the ground state eigenfunction of the harmonic oscillator. The eigenfunctions of the original operator Lj are then given by 2 p˜ j2 x˜j φ˜ mj ,nj = exp − − (30) φmj ,nj . 4 4 Next, we will show that solution of Eq. (20) can be expanded as ρw (t) =

ρ˜m,n

m,n

j

1 φ˜ mj nj . mj !nj !

(31)

By Eq. (20) and multiplying by substituting2 it into 2 ˜ ˜ φ exp( x /4 + p /4) (the left right eigenfunction) m ,n j j j j j from the left and doing the integration, we can obtain the equations of motion for the expansion coefficients: i ∂ ρ˜m,n (t) = − [L, ρ˜m,n ] + [i(mj − nj )ωj ]ρ˜m,n ∂t ¯ j +

i ˆ ρ˜m−j ,n ] mj cj αj(1) [f (q), ¯ j

−i

ˆ ρ˜m−j ,n } mj cj αj(2) {f (q),

j

+

i ˆ ρ˜m+j ,n ] cj αj(1) [f (q), ¯ j

+

1 ˆ ρ˜m,n−j ] nj cj αj(1) [f (q), ¯ j

+

ˆ ρ˜m,n−j } nj cj αj(2) {f (q),

j

−

1 ˆ ρ˜m,n+j ], cj αj(1) [f (q), ¯ j

(32)

where

¯ coth(¯ωj β/2) , 4ωj

1 . 4¯ωj coth(¯ωj β/2) (33) By comparing Eq. (C1) with Eq. (32), we find that the expansion coefficients ρ˜m,n can be calculated using the ADOs ρ m, n αj(1)

=

αj(2)

=

We now present calculations of dissipative quantum dynamics with general spectral densities based on the methods presented in Sec. II, using examples of the spin-Boson dynamics, and the absorption line shape of a model molecular dimer.

A. Spin-Boson dynamics

In the spin-Boson model,44, 45 the system Hamiltonian is Hˆ S = σˆ z + V σˆ x , where σˆ x = |12| + |21| and σˆ z = |11| − |22|. The system operator that couples to the bath is assumed to be qˆs = σˆ z . The initial density matrix of the total system is set to be the factorized form, ˆ ˆ ρT (0) = |11| ⊗ e−β HB /ZB , where ZB = Tre−β HB . 1. Multiple Lorentzian modes

We first consider the case of Ohmic spectral density with exponential cut-off, which has been widely used in previous numerical simulations, where J (ω) = ηωe−ω/ωc .

(35)

A decomposition scheme of J(ω) in Eq. (35) into the form of Eq. (10) has been given by Meier and Tannor,38 where three Lorentzian modes are used. For convenience, the parameters for the three-mode model are given in Table I. Fig. 1 shows the result of the spin-Boson dynamics with the parameters = V = 1, βV = 5.0, η/V = 0.157, and ωc /V = 7.5. Dynamics of this model was first calculated in Ref. 56 using the quasi-adiabatic path integral (QUAPI) method. It has since been used frequently as a numerical benchmark in many studies using both exact and approximate methods. The results in Fig. 1 indicate that, since the systembath coupling is not very strong (the reorganization energy λ ∞ dω = 0.375 V , which is smaller than the cou= π1 0 J (ω) ω pling), the population dynamics converges at the third tier of ADOs, i.e., the sixth order perturbation. In fact, the fourth TABLE I. Parameters for the three Lorentzian mode decomposition of the Ohmic spectral density with exponential cut-off [Eq. (35)] using Eq. (10).38 Mode 1 2 3

pk /(ηωc4 )

k /ωc

k /ωc

12.0677 −19.9762 0.1834

0.2378 0.0888 0.0482

2.2593 5.4377 0.8099


Liu et al.

134106-6

J. Chem. Phys. 140, 134106 (2014)

1 nd

2 order th 4 order th 6 order th 8 order QUAPI

0.4

Error

0.2

0.15 0.1

0

0.05

0

0

-0.005 0

10

5

Super Ohmic Individual Lorentzian modes Sum of Lorentzian modes

2

0.6

0.2

J(ω)ωc /η

P1(t)

0.8

-0.05 0

15

10

5

ω/ωc

Vt

order perturbation (truncation at the second tier) results are already very good. Since this is a low temperature problem, a large number (50–60) of Matsubara terms are needed. In this case, application of the Padé approximation57 can substantially reduce the number of expansion terms in Eq. (8), where converged results can be obtained with a total of 10 exponential terms (six terms come from the poles of the three Lorentzian modes, and the other four terms come from the poles of the Boson factor). Next, we show that the multiple Lorentzian mode scheme can be applied to the super Ohmic spectral density with exponential cut-off: J (ω) =

η ω3 −ω/ωc e . 3! ωc2

(36)

The super Ohmic spectral density has been widely used in several recent studies using the second order quantum master equation in combination with the polaron transformations to calculate the excitation energy transfer dynamics.28, 30 To apply the HEOM formalism to treat the case of super Ohmic spectral density, we first fit it to the analytical form in Eq. (10). The simulated annealing method, as used in Ref. 38, is employed to decompose the spectral density in Eq. (36) with four Lorentzian modes. The fitted parameters are listed in Table II. The comparison between the original

FIG. 2. Numerical fit of the super Ohmic spectral density in Eq. (36) using the sum of four Lorentzian modes.

and fitted spectral densities is shown in Fig. 2, where a very good agreement between the fitted and original spectral densities is found. Two examples for the population dynamics with the super Ohmic spectral density are presented in Fig. 3 with parameters V = 0.5, ωc /V = 2.0, η/V = 2.0, βV = 0.5, /V = 1.0; and in Fig. 4 with parameters V = 0.5, ωc /V = 2.0, η/V = 6.0, βV = 0.5, and /V = 1.0. The dynamics were found to converge with sixth to eighth order perturbations. 2. Combined Debye and oscillatory Debye modes

Examples for the spin-Boson dynamics with combined Debye and oscillatory Debye modes using Eq. (11) are presented in Figs. 5 and 6. The spectral density is taken from Ref. 41, which was calculated numerically from excitation energy fluctuations of bacteriochlorophyll a (Bchla) molecules in the B850 ring of the light harvesting II (LH2) complexes. This spectral density involves two Debye modes and 20 oscillatory Debye modes. The other parameters for the spin-Boson problem are = 0, T = 300 K. Figs. 5 and 6 show results with 1 nd

2 order th 4 order th 6 order th 8 order

0.8

P1(t)

FIG. 1. Evolution of the population on state |1 in the spin-Boson model with the Ohmic spectral density in Eq. (35). The parameters are: = V = 1, βV = 5.0, η = 0.157, and ωc /V = 7.5. The errors in the lower panel are calculated using the converged 12th order result as the reference.

20

15

0.6 0.4

Mode 1 2 3 4

pk /(ηωc4 ) 10.6 7.00 − 0.300 − 6.47

k /ωc 2.48 4.23 0.0581 9.87

k /ωc 2.20 2.25 1.32 3.49

Error

0.2 TABLE II. Parameters for the four Lorentzian mode decomposition of the super Ohmic spectral density with exponential cut-off [Eq. (36)] using Eq. (10), fitted from the simulated annealing method.

0

-0.01 0

5

t

10

FIG. 3. Evolution of the population on state |1 in the spin-Boson model with the super Ohmic spectral density. The parameters are: V = 0.5, ωc/V = 2.0, η/V = 2.0, βV = 0.5, /V = 1.0. The errors in the lower panel are calculated using the converged 12th order result as the reference.


134106-7

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

1

1 nd

P1(t)

0.8 0.6 0.4

nd

2 order th 4 order th 6 order th 8 order th 10 order

0.8

P1(t)


0.6 0.4

0.2

Error

Error

0.05 0

-0.05 -0.1 0

5

10

t

0.01 0

-0.01 0

100

200

300

400

t (fs)

FIG. 4. Same as Fig. 3, except that η/V = 6.0.

FIG. 6. Same as Fig. 5, except that V = 219 cm−1 .

different coupling constants, V = 110 cm−1 and 219 cm−1 , respectively. It can be seen that the results converge with eighth order perturbations.

number of bath modes. But it is found not to be a problem for the simulation time length in Fig. 7, and the result shows that the HEOM using undamped harmonic oscillator modes is a viable simulation method.

3. Discrete harmonic oscillator modes

An example for the HEOM with discrete harmonic oscillator modes is presented in Fig. 7. We use the same parameters as in Fig. 1. The spectral density is decomposed to 60 harmonic oscillator modes, whose frequencies are evenly distributed in the regime between 0 and 4ωc . We note that there are other popular bath discretization methods that have been used to improve the efficiency of the multi-configuration time-dependent Hartree (MCTDH) and the semiclassical simulations of condensed-phase quantum dynamics.58, 59 These methods may also be combined with the discrete harmonic mode decomposition method to further reduce the computational costs. Comparing to the decomposition schemes using dissipative modes as described in Secs. II A and II B, a drawback of the discrete harmonic oscillator method is that recurrence might happen at sufficiently long simulation time with finite

B. Absorption line shape of model dimers

We now calculate the absorption line shape of a model molecular aggregate dimer, to investigate the effects of different types of spectral density. The Frenkel exciton model for two coupled two-level systems is considered. The total Hamiltonian is written as Hˆ tot = Hˆ e + Hˆ ph + Hˆ el−ph .

The first term in the above Eq. (37) denotes the electronic Hamiltonian. When calculating the absorption spectra, only the ground and first excitonic states are considered, Hˆ e = 0 |00| +

2

j |j j | − J (|12| + |21|) ,

1 nd

0.8 0.7 0.6

nd

2 order th 4 order th 6 order QUAPI

0.8 0.6

P1(t)


0.9

P1(t)

(38)

j =1

1

Error

(37)

0.5 0.02

0.4

0.01

0.2

0

-0.01 0

100

200

300

400

t (fs) FIG. 5. Evolution of the population on state |1 in the spin-Boson model with spectral density described in Sec. III A 2. Other parameters are: = 0, V = 110 cm−1 , and T = 300 K. Errors in the lower panel are calculated using the converged 12th order result as the reference.

0 0

10

5

Vt FIG. 7. Evolution of the population on state |1 in the spin-Boson model with the Ohmic spectral density, using the HEOM equation for undamped harmonic oscillators as described in Sec. II C. The parameters are the same as those in Fig. 1.


134106-8

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

where |j represents the state where only the jth molecule is excited; j is the adiabatic transition energy of the jth molecule, and J describes the electronic coupling between the two molecules. The case J > 0 corresponds to a typical J-aggregate.60 The second term Hˆ ph in Eq. (37) is the Hamiltonian of the phonon bath, and is treated as a set of harmonic oscillators. The third term Hˆ el−ph describes the electron-phonon coupling Hˆ el−ph =

2 j =1

cj k xˆk |j j | ,

(39)

k

where the cjk s describe the coupling of the excitation of the jth molecule to the harmonic bath. We further assume that the bath correlation functions for all the molecules have the same spectral density. The absorption line shape can be obtained by calculating the dipole-dipole autocorrelation,17, 52

+∞ 1 I (ω) ∝ dteiωt μ(t) ˆ μ(0) ˆ (40) g, 2π −∞ where ˆ B [e−i H t μρ ˆ g ei H t ]], μ(t) ˆ μ(0) ˆ g = TrS [μTr ˆ

ˆ

(41)

with the total dipole operator μˆ defined as μˆ = 2j =1 μj (|0j | + |j 0|). We have calculated the absorption line shape with the Debye, Ohmic, and super Ohmic spectral densities with exponential cut-off. The Debye spectral density is written as J (ω) =

ηγ ω , + γ2

(42)

ω2

while the Ohmic and super Ohmic spectral densities are described by Eqs. (35) and (36). The three different spectral densities were set to have the same cut-off frequency ωc /J = 1, reorganization energy λ = 2J (the corresponding η parameters are: ηdebye = 4J, ηohmic = 2π J, and ηsohmic = 6π J, respectively), and βωc = 1. It can be seen in Fig. 8 that the line shapes from the Debye and

0.4

Debye Ohmic super-Ohmic

I(ω) (a.u.)

0.3

0.2

0.1

0

-5

0

5

(ω−ε)/J FIG. 8. Calculated absorption spectra of model dimer with different forms of spectral densities. The parameters are: 1 = 2 = , λ/J = 2.0, ωc /J = 2.0, and βωc = 1.0.

Ohmic spectral density are quite similar, while that of the super Ohmic spectral density shows a rather different feature with a strong zero phonon line (ZPL) at the low frequency side of the absorption spectra.

IV. CONCLUSIONS AND DISCUSSIONS

In summary, we have shown that, in combination with different bath spectral density decomposition schemes, as well as some new numerical algorithms like the dynamic filtering and Padé approximation for the bath correlation function, the HEOM method can be efficiently applied to study dissipative quantum dynamics with arbitrary bath spectral densities. The high order perturbation calculations presented in the HEOM formulation are well suited for simulations in the intermediate coupling regime, since convergence can be archived at a not very high order of truncation. In comparison to other high order perturbation methods, the HEOM formalism also has the advantage to assess convergence of the results at certain order of perturbation, and to calculate system-bath correlations through the methods presented in Secs. II D and II E. Within the harmonic bath model, we have also shown that the HEOM formalism is capable to treat undamped harmonic oscillator modes explicitly. In such cases, the HEOM can be regarded as a basis expansion of the quantum Liouville equation expressed in the partial Wigner transformation. Although unlike the case of using dissipative modes, the efficiency of the HEOM formalism for the undamped oscillator modes may or may not be better than the conventional basis expansion methods. In this work, we did not try to directly retrieve the Wigner distribution for the discrete modes, since there are many harmonic modes decomposed from the continuum bath model. The new approach could be very useful to analyze system-bath correlations in quantum dynamics in cases where only a few number of undamped harmonic oscillator modes are involved. Finally, we have a comment on the spectral density decomposition methods in Eq. (10). Although various types of spectral density can be fitted using the different schemes, as shown in previous study38 and in the paper, there are cases where the asymptotic properties of the spectral density or correlation function are not preserved by the bath decomposition schemes in Eqs. (10) and (11). For example, for the sub-Ohmic spectral density,61, 62 the reorganization energy diverges. However, the expansion in either Eq. (10) or Eq. (11) does not have this property. So caution should be used when trying to investigate dynamic behaviors related to such asymptotic properties when using the fitted spectral densities, where new decomposition schemes may be needed.

ACKNOWLEDGMENTS

This work is supported by NSFC (Grant Nos. 91027015 and 21290194) and the 973 program (Grant Nos. 2011CB808502 and 2013CB933501).


Liu et al.

134106-9

J. Chem. Phys. 140, 134106 (2014)

APPENDIX A: SUM OF LORENTZIAN MODES

Using Eq. (10), the correlation function in Eq. (8) can be written into the following form,37, 38 C(t) =

N

[dj + e−( j +ij )t + dj − e−( j −ij )t ] +

j =1

∞

dk e−νk t ,

N 2pj νk 2i . dk = − J (−iνk ) = − 2 2 2 β β j + j − νk2 + 4νk2 2j j =1 (A3)

k=1

(A1) where

dj ±

ν k = 2π k/β are the Matsubara frequencies, and the corresponding coefficients are

β(j ∓ i j ) pj coth ±1 , = 8j j 2

(A2)

We note that, instead of using the Matsubara expansion, it is also possible to use the Padé approximation for expanding the Boson factor in Eq. (8),57, 63, 64 which is found to be very efficient at low temperatures.57 Details of this scheme can be found in Ref. 57, and will not be presented in this work. Accordingly, the HEOM can be expressed as

⎡ ⎤ N n ∞ ∂ m,n i m,n − ⎣ (mj + nj ) j + i (mj − nj )j + jk γk ⎦ ρJm,n ρ (t) = − HS , ρJ ∂t J ¯ j =1 j =1 k=1

i i i m− m,n− j ,n j ˆ ρJm,n m j + − nj − − dk jk f (q), − j ρJ j ρJ k ¯ j =1 ¯ j =1 ¯ k=1

−

N ∞ i i m+ ,n m,n+ ˆ ρJ j + ρJ j − ˆ ρJm,n f (q), f (q), , + k ¯ ¯ j =1 k=1

N

N

where m, n denote the primary pseudomodes associated with dj± , J denotes the Matsubara terms, and ± j ρ

∞

−

pj β(j ∓ i j ) ˆ ρ] [f (q), = coth 8j j 2 pj {f (q), ˆ ρ} . ± 8j j

given by C(t) =

Na

bj e−γa,j t +

j =1

+

(A5)

The above Eq. (A4) is just an extension of the previous single mode Lorentzian result65 to multiple modes. We further note that the current formalism is also different from the equations originally derived for the second order perturbation GQME by Meier and Tannor,38 as well as the high order perturbation form obtained by Kleinekathöfer and coworkers,37 where the HEOM is derived using the stochastic method.32 The difference is that, in the current formulation, the real and imaginary parts of the bath correlation function with the same exponential constant are treated together rather than separately. Although simulations using the two different forms of HEOM should give the same results, the form in Eq. (A4) needs only half number of modes to represent the bath correlation functions, which is an advantage in numerical simulations.

APPENDIX B: SUM OF DEBYE AND OSCILLATORY DEBYE MODES

With the bath decomposition using Eq. (11), the quantum bath correlation functions expanded in the form of Eq. (8) are

(A4)

∞

Nb

dj + e−(γb,j +iωb,j )t + dj − e−(γb,j −iωb,j )t

j =1

dk e−νk t ,

(B1)

k=1

where bj =

dj +

dj −

β¯γa,j ηa,j γa,j cot −i , 2 2

(B2a)

β¯(ωb,j − iγb,j ) ηb,j (ωb,j − iγ ) coth −i , = 4 2 (B2b)

β¯(ωb,j + iγb,j ) ηb,j (ωb,j + iγ ) coth − i , (B2c) = 4 2

dk =

2ηa,j γa,j νk 2 2 β νk − γa,j 2 2 ν k ωb,j − νk2 + γb,j 4ηb,j γb,j − . 2 2 2 2 β ωb,j − νk2 + γb,j + 4νk2 ωb,j

(B2d)


134106-10

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

The corresponding HEOM can be obtained ⎡ Na Nb ∂ m,n i m,n −⎣ rj γa,j + (mj + nj )ωb,j ρr,J (t) =− Hs , ρr,J ∂t ¯ j =1 j =1 − i(mj − nj )γb,j +

∞

m,n jk γk ρr,J

a a

i i ˆ ρrm,n f (q), rj j ρrm,n − − + ,J ,J j j ¯ j =1 ¯ j =1

J

J

−

i ¯

j =1 ∞

−

m,n− j

nj j − ρr,J

−

i ¯

Jb

∂ ρm,n (t) ∂t

−i

Jb Jb

i i m− ,n m+ ,n ˆ ρr,Jj − f (q), mj j + ρr,Jj − ¯ j =1 ¯ j =1 Jb

The HEOM with the discrete harmonic oscillator modes in the spectral density is given by

N i = − [L, ρm,n ] + [i(mj − nj )ωj ]ρm,n ¯ j =1

k=1

−

APPENDIX C: USING UNDAMPED HARMONIC OSCILLATOR MODES

m,n+ ˆ ρr,J j f (q),

j =1

j =1

+i

−i

∞

−i

where the subscripts r are used for the Debye modes, m and n are for the under damped Debye modes, and J are for the Matsubara modes. The other labels with a ± sign differ from {r, m, n, J} denote changes of the specified index to one-level higher or lower. The operators j , j + , and j − are defined as

j + ρ =

j − ρ =

ηa,j γa,j 2

βγa,j ˆ ρ] − i{f (q), ˆ ρ} , cot [f (q), 2 (B4a)

ηb,j (ωb,j − iγ ) 4 β¯(ωb,j − iγb,j ) ˆ ρ] − i{f (q), ˆ ρ} , × coth [f (q), 2 (B4b)

ηb,j (ωb,j + iγ ) 4 β¯(ωb,j + iγb,j ) ˆ ρ] − i{f (q), ˆ ρ} . [f (q), × coth 2 (B4c)

We also note that although Eqs. (B2a)–(B2d) employ the Matsubara expansion scheme, the Padé scheme to decompose the Boson factor57 can also be applied to increase the efficiency at low temperatures.

4ωj

N nj cj2 j =1

i i m,n m,n ˆ ρr,J ˆ ρr,J f (q), dk jk f (q), − − + , k k ¯ k=1 ¯ k=1

4ωj

N mj cj2 j =1

4ωj

N nj cj2 j =1

(B3)

j ρ =

N mj cj2

4ωj

coth

β¯ωj ˆ ρm−j ,n ] [f (q), 2 i ˆ ρm+j ,n ] [f (q), ¯ j =1 N

ˆ ρm−j ,n }− {f (q),

coth

β¯ωj ˆ ρm,n−j ] [f (q), 2 i ˆ ρm,n+j ], (C1) [f (q), ¯ j =1 N

ˆ ρm,n−j } − {f (q),

where the subscripts m and n denote the modes with exponential coefficients {iωj } and {−iωj }. APPENDIX D: DERIVATION OF EQ. (18)

In this Appendix, we present the detailed derivation of the expansion coefficients aj in Eq. (18) regarding the distribution function of the collective bath coordinate. In Ref. 46, it has been shown that the expectation value of powers concerning the collective bath coordinate can be calculated using the ADOs: (n+1) j! ρn , Cj X(n+1) = k nk ! j ≥0 j th order (D1) nk = j, n + 1 ≥ j, where (n−1) Cj(n+1) = −Cj(n) , −1 + nA0 Cj

(D2)

with C00 = 1. In the above Eqs. (D1) and (D2), j = 1, 3, . . . , n + 1 when n is even, and j = 0, 2, 4, . . . , n + 1 when n is odd. On the other hand, according to Eq. (15), the nth order expectation value of X can be written as Xn = aj Kjn , (D3) j

where

Kjn =

dXXn j (X),

(D4)

with j (X) defined in Eq. (16). Using the technique of integration by parts, and the recurrence relation for the Hermite polynomial Hj (x) = 2j Hj −1 (x), one gets n−2 Kjn = j A0 Kjn−1 . (D5) −1 + (n − 1)A0 Kj


134106-11

Liu et al.

J. Chem. Phys. 140, 134106 (2014)

It can be shown that recurrence formulas, Eqs. √ √ (D2) and (D5), are equivalent if we set Kjn = (− A0 )j j !Cjn . Hence, by comparing Eqs. (D1) and (D3), we can obtain the expansion coefficients as (−1)j aj = j !A0 j

j th order

j! ρn , k nk !

nk = j,

(D6)

which is just Eq. (18). 1 V.

May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed. (Wiley-VCH, Weinheim, 2011). 2 P. F. Barbara, T. J. Meyer, and M. A. Ratner, J. Phys. Chem. 100, 13148 (1996). 3 Y.-C. Cheng and G. R. Fleming, Annu. Rev. Phys. Chem. 60, 241 (2009). 4 D. Beljonne, C. Curutchet, G. D. Scholes, and R. J. Silbey, J. Phys. Chem. B 113, 6583 (2009). 5 T. Renger, V. May, and O. Kühn, Phys. Rep. 343, 137 (2001). 6 D. Abramavicius, B. Palmieri, D. V. Voronine, F. Sanda, and S. Mukamel, Chem. Rev. 109, 2350 (2009). 7 V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Oliver, R. Silbey, and J. L. Brédas, Chem. Rev. 107, 926 (2007). 8 R. A. Marcus and N. Sustin, Biochim. Biophys. Acta 811, 265 (1985). 9 R. A. Marcus, Rev. Mod. Phys. 65, 599 (1993). 10 T. Förster, Ann. Physik 437, 55 (1948). 11 A. G. Redfield, IBM J. Res. 1, 19 (1957). 12 W. T. Pollard, A. K. Felts, and R. A. Friesner, Adv. Chem. Phys. 93, 77 (1996). 13 Y. J. Yan and R. X. Xu, Annu. Rev. Phys. Chem. 56, 187 (2005). 14 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234110 (2009). 15 A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009). 16 M. Yang and G. R. Fleming, Chem. Phys. 282, 163 (2002). 17 L.-P. Chen, R.-H. Zheng, Q. Shi, and Y.-J. Yan, J. Chem. Phys. 131, 094502 (2009). 18 L.-P. Chen, R.-H. Zheng, Y.-Y. Jing, and Q. Shi, J. Chem. Phys. 134, 194508 (2011). 19 B. B. Laird, J. Budimir, and J. L. Skinner, J. Chem. Phys. 94, 4391 (1991). 20 D. R. Reichman and R. J. Silbey, J. Chem. Phys. 104, 1506 (1996). 21 S. Jang, J. Cao, and R. J. Silbey, J. Chem. Phys. 116, 2705 (2002). 22 L.-P. Chen and Q. Shi, J. Chem. Phys. 130, 134505 (2009). 23 D. Wang, L. P. Chen, R. H. Zheng, L. J. Wang, and Q. Shi, J. Chem. Phys. 132, 081101 (2010). 24 L.-P. Chen, R.-H. Zheng, Q. Shi, and Y.-J. Yan, J. Chem. Phys. 132, 024505 (2010). 25 M. Sparpaglione and S. Mukamel, J. Chem. Phys. 88, 3263 (1988). 26 A. A. Golosov and D. R. Reichman, J. Chem. Phys. 115, 9848 (2001). 27 A. Pomyalov and D. J. Tannor, J. Chem. Phys. 123, 204111 (2005).

28 S. Jang, Y.-C. Cheng, D. R. Reichman, and J. D. Eaves, J. Chem. Phys. 129,

101104 (2008). Kolli, A. Nazir, and A. Olaya-Castro, J. Chem. Phys. 135, 154112 (2011). 30 H.-T. Chang and Y.-C. Cheng, J. Chem. Phys. 137, 165103 (2012). 31 Y. Tanimura and R. K. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989). 32 Y. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395, 216 (2004). 33 R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo, and Y.-J. Yan, J. Chem. Phys. 122, 041103 (2005). 34 A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn. 74, 3131 (2005). 35 Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). 36 Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130, 084105 (2009). 37 M. Schröder, M. Schreiber, and U. Kleinekathöfer, J. Chem. Phys. 126, 114102 (2007). 38 C. Meier and D. Tannor, J. Chem. Phys. 111, 3365 (1999). 39 U. Kleinekathöfer, J. Chem. Phys. 121, 2505 (2004). 40 M. Schröder, U. Kleinekathöfer, and M. Schreiber, J. Chem. Phys. 124, 084903 (2006). 41 C. Olbrich and U. Kleinekathöfer, J. Phys. Chem. B 114, 12427 (2010). 42 C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett. 3, 2828 (2012). 43 A. O. Caldeira and A. J. Leggett, Physica A 121, 587 (1983). 44 A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987). 45 U. Weiss, Quantum Dissipative Systems, 2nd ed. (World Scientific, London, 1999). 46 L.-L. Zhu, H. Liu, and Q. Shi, J. Chem. Phys. 137, 194106 (2012). 47 E. Wigner, Phys. Rev. 40, 749 (1932). 48 J. Strümpfer and K. Schulten, J. Chem. Phys. 131, 225101 (2009). 49 T. Joo, Y. Jia, J. Yu, M. J. Lang, and G. R. Fleming, J. Chem. Phys. 104, 6089 (1996). 50 M. Yang, R. Agarwal, and G. R. Fleming, J. Photochem. Photobiol. A 142, 107 (2001). 51 A. Ishizaki and Y. Tanimura, Chem. Phys. 347, 185 (2008). 52 S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford, New York, 1995). 53 R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919 (1999). 54 K. Ando and M. Santer, J. Chem. Phys. 118, 10399 (2003). 55 Q. Shi and E. Geva, J. Chem. Phys. 121, 3393 (2004). 56 N. Makri and D. Makarov, J. Chem. Phys. 102, 4600 (1995). 57 J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134, 244106 (2011). 58 H. Wang, J. Chem. Phys. 113, 9948 (2000). 59 H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115, 2979 (2001). 60 J-Aggregates, edited by T. Kobayashi (World Scientific, Singapore, 1996). 61 R. Bulla, N. Tong, and M. Vojta, Phys. Rev. Lett. 91, 170601 (2003). 62 H. Wang and M. Thoss, Chem. Phys. 370, 78 (2010). 63 T. Ozaki, Phys. Rev. B 75, 035123 (2007). 64 J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133, 101106 (2010). 65 M. Tanaka and Y. Tanimura, J. Chem. Phys. 132, 214502 (2010). 29 A.


Analytic representations of bath correlation functions for ohmic and superohmic spectral densities using simple poles.

Dissipative dynamics of a driven quantum spin coupled to a bath of ultracold fermions.

Markovian master equation for a classical particle coupled with arbitrary strength to a harmonic bath.

Bed bath.

Non-adiabatic holonomic quantum computation in linear system-bath coupling.

Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities.

Water-bath temperature controller.

Synthetic cathinones ("bath salts").

Sauna bath burn.

Bath-taking made easy.

Noncanonical statistics of a finite quantum system with non-negligible system-bath coupling.

Hostel for Girls at Bath.

Mineral bath therapy in arthritis.

Study of the morphology of ZnS thin films deposited on different substrates via chemical bath deposition.

Transcriptomic profiling of Methylococcus capsulatus (Bath) during growth with two different methane monooxygenases.

Evaluation of a new denture bath solution.

Suppression of nuclear spin bath fluctuations in self-assembled quantum dots induced by inhomogeneous strain.

Bath-induced paroxysmal events - Notes on a case report.

The Great Arcanum:-The Turkish Bath.

Medicines advertised in eighteenth-century Bath newspapers.

Bubble bath burns: an unusual case.

Unexplained acute fever after a hot bath.

Tissue bath improvements for the Oxford Vibratome.

Bath salts and other emerging toxins.