Article pubs.acs.org/JPCA
Chebyshev Expansion Applied to Dissipative Quantum Systems Bogdan Popescu,†,‡ Hasan Rahman,† and Ulrich Kleinekathöfer*,† †
Department of Physics and Earth Sciences, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany ABSTRACT: To determine the dynamics of a molecular aggregate under the influence of a strongly time-dependent perturbation within a dissipative environment is still, in general, a challenge. The time-dependent perturbation might be, for example, due to external fields or explicitly treated fluctuations within the environment. Methods to calculate the dynamics in these cases do exist though some of these approaches assume that the corresponding correlation functions can be written as a weighted sum of exponentials. One such theory is the hierarchical equations of motion approach. If the environment, however, is described by a complex spectral density or if its temperature is low, these approaches become very inefficient. Therefore, we propose a scheme based on a Chebyshev decomposition of the bath correlation functions and detail the respective quantum master equations within second-order perturbation theory in the environmental coupling. Similar approaches have recently been proposed for systems coupled to Fermionic reservoirs. The proposed scheme is tested for a simple two-level system and compared to existing results. Furthermore, the advantages and disadvantages of the present Chebyshev approach are discussed.
I. INTRODUCTION To tackle the external influence of, e.g., a protein environment on the excitation energy transfer in systems, such as pigment− protein complexes, various methods based on the reduced density operator (RDO) are available.1−3 The time evolution of the latter can, for instance, be determined using quantum master equations (QMEs).1−3 However, these have severe limitations as they are commonly derived in the lowest-order perturbation theory. To overcome this bottleneck, the hierarchical equations of motion (HEOM) approach is one possible and popular possibility.4−8 The HEOM theory was initially formulated for dissipative systems by Tanimura and Kubo4 and extended to Fermionic systems by Yan and coworkers.5 By now it has already been applied to a wide range of bosonic6,9,10 and Fermionic scenarios.5,11−13 Most notably, for the Fermionic case, its second-tier truncation is exact for noninteracting particles and was connected to Keldysh nonequilibrium Green’s functions, whereas the first-tier corresponds to time-nonlocal master equations.13 However, HEOM calculations for large systems are still challenging, merely due to their large computational demands.7,8 Some RDO theories such as the HEOM approach involve a so-called multi-Lorentzian parametrization of the reservoir spectral densities14 as well as a Matsubara or equivalent expansion of the Fermi function.15,16 These expansions turn the bath correlation functions (BCFs) or, equivalently, the self-energies of the leads into a sum of weighted exponentials. Although feasible in principle, the reconstruction of environmental spectral densities in terms of Lorentzian functions can be rather cumbersome in practice.17 For example, as discussed in ref 18 for the case of graphene nanoribbons, the Lorentzian fitting of a given lead self-energy can be pursued by employing © 2016 American Chemical Society
a least-squares solver with the goal of minimizing deviation functions between the fitted object and its Lorentzian reconstruction. Moreover, in particular cases such as lightabsorbing aggregates, where the spectral densities show a nonlinear behavior for small frequencies, multi-Lorentzian fitting fails in that range. Therefore, new types of parametrizations were proposed19 using other classes of fit functions better suited to mimic superohmic behavior. The fitting procedure becomes even more complicated when we try to embed experimental or numerically obtained spectral densities, which, at least for certain systems, are known to contain a strongly varying peak structure.20,21 For Fermionic reservoirs, Tian and Chen11,12 employed the Jacobi−Anger identity to expand the HEOM (truncated at the second-tier) in terms of Chebyshev polynomials. This scheme has the advantage that the correlation function does not have to be a weighted sum of exponentials, which also implies that very low temperatures can be used without any additional problems. At the same time, Chebyshev expansions do have the drawback that the simulation time is limited by the truncation order of the series that is infinite in principle. In ref 11 the auxiliary density operators (ADOs) and the lead spectral density are decomposed in terms of Chebyshev polynomials. Recently, we took up the idea of using a Chebyshev expansion but used it directly to decompose the BCF, which is much more in the spirit of the original perturbative or HEOM schemes.22−24 To this end, a scheme for molecular wire scenarios, i.e., for Special Issue: Ronnie Kosloff Festschrift Received: December 14, 2015 Revised: February 4, 2016 Published: February 4, 2016 3270
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A
reservoir coupling Hamiltonian HSR can be written as a sum of products of system and bath operators. The latter are denoted by Φj, whereas the former are supposed to be site diagonal for excitonic systems3
Fermionic leads, has been detailed to second-order in the wirelead coupling.25 Compared to the scheme by Tian and Chen,11,12 the scheme by the present authors only includes a single expansion instead of a 2-fold one. So there is only one truncation step that might turn out to be advantageous. The goal of the present contribution is to adapt the QME derived in ref 25 to the case of a system in contact with bosonic reservoirs. Note that polynomial expansions and especially Chebyshev expansions have been applied extensively also in the context of wave packet propagation26−29 also for inhomogeneous30 or explicitly time-dependent cases.31 For open quantum systems Chebyshev polynomials32 as well as their generalizations, the Faber polynomials have been employed, as have been Newton polynomials.33 Recently, a decomposition scheme termed timeevolving density with orthogonal polynomials algorithm (TEDOPA) has been reported which is able to treat dissipative systems with arbitrary spectral densities.34 The article is organized as follows. In the subsequent section the tight-binding modeling, typically used for open quantum systems is discussed whereas in section III the Chebyshev expansion is detailed. In section IV we employ the Chebyshev expansion to derive QMEs in second-order perturbation theory. Finally, in section V, numerical results are shown and discussed. In particular, we apply the new equations to a typical benchmarking two-level model and compare the results with other standard methods such as the HEOM.
HSR =
j
1j(ω) =
i≠j
(3)
ξ
cjξ 2 2mξ ωξ
δ(ω − ωξ) (4)
For a set of bath oscillators the spectral density consists of a sequence of δ-functions, assumed sufficiently dense to form a continuous spectrum. Quantitatively, 1j(ω) is related to the bath correlation functions by3 ℏ π ℏ = π ℏ = π
Cj(t ) =
∞
∫−∞
dω
e−iωt 1j(ω) 1 − e−βℏω
∞
∫−∞ dω e−iωt 1j(ω)[1 + n(ℏω)] ∫0
∞
dω 1j(ω)[cos(ωt ) coth(β ℏω /2) − i sin(ωt )] (5)
where β = 1/kBT denotes the inverse temperature and n(ℏω) the Bose−Einstein distribution.
III. CHEBYSHEV DECOMPOSITION OF CORRELATION FUNCTIONS In this section, we sketch the line of reasoning that we followed in ref 25 to expand the BCF in terms of Chebyshev polynomials, however, adapted here to bosonic environments. To make use of this type of orthogonal polynomial, a mapping of the spectral range belonging to the Hamiltonian to their interval of definition is in order, i.e., to the interval [−1, +1]. To this end, one first needs to restrict the infinite integration range in eq 5 to a finite range [ωmin, ωmax]. Then, by introducing the dimensionless variable x = (ω − ω̅ ) /Ω, where ω̅ = (ωmax + ωmin)/2 and Ω = (ωmax − ωmin)/2, the integral expression for the correlation function turns into
(1)
Note that in the case of exciton transfer, which is often described by the above Hamiltonian, the exciton energies can be determined from electronic structure calculations together with the electronic couplings and are mapped onto the above tight-binding model.20 The bosonic reservoir is treated as an infinite collection of displaced harmonic oscillators, linearly coupled to the system degrees of freedom and acting as a heat bath ⎛ p2 mξ ωξ 2xξ 2 ⎞ ξ ⎜ ⎟ HR = ∑ ⎜ + ⎟ 2mξ 2 ⎠ ξ ⎝
∑ ξ
∑ Ei|i⟩⟨i| + ∑ Vij|i⟩⟨j| i
j
Here the cjξ define the coupling strengths of the bath oscillators to the system states and the system part of the system− reservoir coupling is given by a projector onto site j, i.e., Kj = | j⟩⟨j|. In this form an interaction between the bath and the system at site j exists only if an excitation is present. Also note that we restricted ourselves to situations in which each site has its own thermal bath, which is assumed to be uncorrelated to the bath modes at other sites. The effect of the environment on the system is typically given in form of the spectral density, which contains information concerning the spectrum of the environment and the system−bath coupling
II. MODEL HAMILTONIAN In theories for open quantum systems, the Hamiltonian of a composite system is typically decomposed into three terms according to H = HS + HR + HSR. Herein, the Hamiltonian of the quantum system is denoted by HS(t), the one of the bosonic reservoir, playing the role of the environment, by HR, whereas HSR stands for the system−reservoir coupling Hamiltonian. Because a particular focus is placed on the transfer of a single electronic excitation, we restrict ourselves to the single exciton manifold. Thus, the system basis states |i⟩ assume one exciton at site i, corresponding to a molecular orbital, whereas the other sites are supposed to be in their ground states. The ground state without any excitations |0⟩ is not explicitly considered here. Moreover, the electronic system is modeled by N sites with site energies Ei at site i and electronic couplings Vij between any site i and j, leading to the following tight-binding description for HS HS =
∑ ΦjKj = ∑ ∑ cjξxξ|j⟩⟨j|
Cj(t ) =
ℏΩ π
1
∫−1
dx
e−i(Ωx + ω̅ )t 1j(Ωx + ω̅ ) 1 − e−βℏ(Ωx + ω̅ )
(6)
Next, we expand the Fourier kernel in eq 6 by means of the Jacobi−Anger identity35 given by ∞
(2)
e−iΩxt = J0 (Ωt ) +
In this equation pξ, xξ, ωξ, and mξ denote the momenta, displacements around the equilibrium, frequencies, and masses of the bath oscillators, respectively. Finally, the system−
∑ 2(−i)n Jn(Ωt ) Tn(x) n=1
∀t ∈ , ∀x ∈ [−1, +1] 3271
(7) DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A d ρ ̃ (t ) dt S
with the Chebyshev polynomials Tn and the Bessel functions of the first kind Jn. Subsequently, the BCF can be written as ∞
Cj(t ) =
=−
−iω t
∑ Jn (Ωt )e
1 ℏ2
∫0
t
dτ TrR {[H̃ SR (t ), [H̃ SR (τ ), ρS̃ (t ) ρReq ]]}
n=0
(11)
Cn(t )
×
ρeq R
where denotes the equilibrium density operator of the reservoir, ρS(t) is the system reduced density operator, and the tilde on top of the operator symbols marks the interaction picture. Equation 11 can be expressed in a more convenient form by inserting the expression of the coupling Hamiltonian, eq 3, and subsequently changing to the Schrödinger picture
1j(Ωx + ω ) ℏΩ 1 dx (2 − δ0, n)(− i)n Tn(x) −1 1 − e−βℏ(Ωx + ω) π
∫
Ij , n ∞
=
∑ Cn(t )Ij ,n
(8)
n=0
d i ρ (t ) = − 3S(t ) ρS (t ) dt S ℏ 1 − ∑ [Kj , Λj(t ) ρS (t ) − ρS (t ) Λ†j (t )] ℏ j
where δ0,n denotes the Kronecker delta. In eq 8 the timedependent quantities Cn(t) as well as the time-independent integrals Ij,n were defined. It is crucial to note that the quantities Ij,n need to be computed only once at runtime, whereas the time dependence of the BCF is encapsulated in the simpler expressions Cn(t). Moreover, up to this point there were no assumptions on the form of spectral densities nor on the Bose− Einstein distribution, hence also on the involved temperature. This suggests that such an expansion might be advantageous for embedding complex spectral densities and low temperatures. This fact is substantially in contrast with the exponential decomposition of the BCF, which becomes very demanding for complex spectral densities and/or low temperatures. We would like to note that the numerical evaluation of the integrals Ij,n by standard techniques turns out to be nontrivial due to the strongly oscillating Chebyshev polynomials for larger values of n. However, this can be circumvented by the substitution x = cos(Θ) and the subsequent use of the property Tn(cos(Θ)) = cos(Θ), which renders the integrands to the general form f(x)· cos(ωx), thus amenable to specifically designed integration subroutines, e.g., from the QUADPACK library.36 The key point in deriving a QME in the exponential decomposition of the correlation function22 combined with a complex-pole expansion of the Bose−Einstein distribution is to be able to write the BCF as a weighted sum over exponentials. This decomposition guarantees that time derivatives of partial BCFs can be expressed in terms of themselves. However, to achieve the latter goal, one can instead benefit from the Bessel functions of the first kind entering the time-dependent part of the BCF in eq 8 because these have a similar property35 d Ω Jn (Ωt ) = (Jn − 1(Ωt ) − Jn + 1(Ωt )) dt 2
Consequently, the time derivative of Cn(t) = given by
Here 3S = [HS(t ), •] denotes the Liouville operator and the Λj(t) are auxiliary density operators (ADOs) defined as22 Λ j( t ) =
∫0
t
dτ Cj(t −τ ) US(t ,τ )Kj US†(t ,τ )
∞
Λ j(t ) =
∑ Ij , n ∫ n=0
t
dτ Cn(t −τ ) US(t ,τ )Kj US†(t ,τ ) 0 Λj , n(t )
∞
=
(9)
d d Cn(t ) = −iω̅ Jn (Ωt )e−iω̅ t + e−iω̅ t Jn (Ωt ) dt dt Ω = −iω̅ Cn(t ) + (Cn − 1(t ) − Cn + 1(t )) 2
(13)
where US(t,τ) stands for the time-evolution operator of the system. The QME in eq 11 is given in the so-called time-local (TL) form22 rather than in its time-nonlocal form.14 In ref 22 this TL form has been combined with the exponential decomposition of the BCF leading to the exponential decomposition time-local (EDTL) formalism. Alternatively, we now make use of the Chebyshev analysis of the BCF detailed in the previous section. This then leads to the Chebyshev decomposition time-local (CDTL) formalism, a notation introduced previously.25 When the respective decomposition is converged, the EDTL and CDTL approaches yield the same results as expected, because they are based on the same time-local QME. Using the Chebyshev expansion turns eq 13 into
∑ Ij ,n Λj ,n(t ) n=0
Jn(Ωt)e−iω̅ t
(12)
(14)
This result consists of a time-dependent part of the ADO encapsulated in the new quantities Λj,n(t), which we term partial ADOs. For later usage, one can derive simple master equations for these partial ADOs by directly differentiating their definition expression, eq 14, and by using the relation, eq 10. In addition, the infinite sum over the partial ADOs is truncated at a finite order NC. As detailed in ref 25 and similar to the method of ref 11, truncating at a certain NC limits the approach to a NC-dependent simulation time. An estimate for the number of terms needed to obtain converged results can be given by25,37,38
can be
(10)
which connects the derivative of Cn(t) to itself and its nearestrank neighbors. This expression is the key point in deriving a QME in the next section.
N = ΩT + 10 ln(ΩT )
(15)
If this relation is fulfilled, the uniform Chebyshev approximation shows an exponential convergence leading to very accurate results.27 Summarizing the final set of QMEs for the quantities Λj,n(t) together with the modified forms of the
IV. QUANTUM MASTER EQUATIONS As discussed in previous studies,2,3,22 the second-order timelocal quantum master equation can be written as 3272
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A
Figure 1. Time evolution of the population on the first site that is initially excited for different reorganization energies. The site energies are assumed the same, i.e., ΔE = E2 − E1 = 0.
QMEs for the boundary cases, the zeroth and NCth order, one obtains ⎧ i ⎪ Kj − 3S(t )Λj ,0(t ) − iω̅ Λj ,0(t ) ℏ ⎪ ⎪ − ΩΛj ,1(t ) if n = 0 ⎪ i ⎪− 3 (t )Λ (t ) − iω Λ (t ) ̅ j,n j,n ⎪ ℏ S d Λj , n(t ) = ⎨ Ω dt ⎪ + [Λj , n − 1(t ) − Λj , n + 1(t )] if 0 > n > N C 2 ⎪ ⎪ i ⎪− 3S(t )Λj , NC(t ) − iω̅ Λj , NC(t ) ⎪ ℏ Ω ⎪ + Λj , N − 1(t ) C if n = NC ⎩ 2
numerically or experimentally obtained spectral densities as well as low temperatures.
V. HIERARCHICAL EQUATIONS OF MOTION AND ENSEMBLE-AVERAGED WAVE PACKET DYNAMICS The HEOM scheme can only be employed in cases where the correlation function can be written as a weighted sum of exponentials, implying that it uses the exponential decomposition.7,10,24 The EDTL scheme22 is actually equivalent to the first tier HEOM formalism if properly terminated. This restriction implies that the spectral density can be written as a sum of Lorentzian functions and that the Fermi function is written as a series expansion such as the Matsubara series. To this end, auxiliary density operators take into account the nonMarkovian dynamics induced by the environment. The exponential decomposition together with the hierarchical scheme leads, in principle, to an infinite number of terms that can be truncated at a finite level, however, depending on the required accuracy. Appropriate approximations can be employed in the case of low temperatures.39,40 In the present case a code described earlier7 has been used to obtain results converged with respect to the number of auxiliary terms. The numerical findings coincide with results already detailed in refs 41, 42, and 43 where also some more details on the convergence have been reported. The other method to which we compare the present CDTL results is the ensemble-averaged wave packet dynamics in its classical-path based version. In this approach, other than in the standard Ehrenfest dynamics as, for example, employed in ref 44, the effect of the system on the bath is neglected. In other words, the bath always stays in equilibrium if it is initially in equilibrium. This assumption is common in perturbative density matrix schemes.3 At the same time, this assumption clearly indicates that the approach will not be able to handle strongly coupled system−bath configurations. In the ensembleaveraged wave packet dynamics, the effect of the environment is
(16)
The above system along with the summation rule (14) and the equation of motion for the RDO, eq 12, form a complete set of coupled differential equations, amenable to numerical treatment. A note is in order concerning the Markovian limit of the obtained set of QMEs. In the present time-local approach, the Markovian limit is identical to moving the upper integration limit in eq 13 or 14 to infinity, i.e., using Λj(t=∞) while propagating the density matrix ρS. Because the time evolution of the Λj(t) is independent of ρS, eq 16 can be used to first determine Λj(t=∞) and then be employed in eq 12 to propagate ρS. The non-Markovian results together with this Markovian limit are tested below for a simple model, i.e., a twolevel dissipative system, showing its performance, in comparison to the well-established HEOM theory. It is worth stressing again the advantage of the present method over previous, similar approaches,23 i.e., the fact that the Bose−Einstein distribution and the bath spectral densities only enter the time-independent quantities Ij,n, which need to be computed only once, namely at the beginning of the simulation. Therefore, the present scheme has no restriction, neither on the form of the respective spectral densitiy nor on the system temperature, thus being able to cope with 3273
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A
Figure 2. Same as in Figure 1 but for unequal site energies ΔE = E2 − E1 = 100 cm−1.
simulation parameters were chosen to be V = 100 cm−1 for the electronic coupling, T = 300 K for the temperature, and 1/γi = 100 fs for the correlation time. Moreover, we distinguish two cases; i.e., in Figure 1 equal site energies E1 = E2 are considered, whereas Figure 2 shows the case of unequal site energies, i.e., an energy splitting of ΔE = E2 − E1 = 100 cm−1. For both cases, i.e., Figures 1 and 2, the reorganization energy, λi, assumes four different values, indicated in each panel spanning the system−bath coupling regime from low to high values. The plots show the population dynamics of the first site, i.e., element (1, 1) of the density operator, given by the present Chebyshev QME as compared to its counterparts computed by the HEOM scheme. The HEOM results including details of their calculation have been reported in ref 43 already. The important point for the present study is that these findings can serve as benchmark results because they are converged in the system−bath coupling. In Figure 1 the results for the case of equal site energies are shown. For the lowest reorganization energy (i.e., λ = 2 cm−1) the agreement with the HEOM results is good, which points out that the present scheme correctly describes the underlying dynamics. Moving on to intermediate values of λ, i.e., Figure 1b,c, the damped oscillations in the Chebyshev dynamics start to differ from those of the HEOM results mainly in amplitude. Moreover, a slight shift in frequency can be observed as well. This finding indicates that this medium magnitude range of the system−bath coupling already marks the start of the breakdown of second-order perturbation theory, which is the basis of the CDTL QME. However, it is interesting to note that the scheme still correctly describes the long-time limit that corresponds to equal population of both sites for equal site energies. Using the CDTL approach, the correct thermal state is also obtained for the largest reorganization energy of λ = 500 cm−1 in Figure 1. This second-order scheme reaches this limit, however, much too fast compared to the converged HEOM results. The case of unequal site energies is shown in Figure 2. The agreement between the two approaches is quite similar to that for the case
given in terms of fluctuations in the site energies ΔEm(t), which are determined as specific matrix elements of bath operators. As a starting point, one determines time-dependent trajectories of site energies including the ΔEm(t) and couplings. Subsequently, the time-dependent Schrödinger equation is solved for the state vector |ΨSα(t)⟩ of ensemble member α along N pieces of the trajectory or an ensemble of trajectories. Finally, an ensemble average over the different pieces of the trajectory has to be performed to obtain the expectation value of an observable A leading to ⟨A(t )⟩ =
1 Nα
Nα
∑ ⟨ΨSα(t )|A|ΨSα(t )⟩ = ⟨ΨS(t )|A|ΨS(t )⟩ α=1
(17)
A comparison of the HEOM and ensemble-averaged wave packet has been reported earlier,43 though not with a secondorder perturbative treatment.
VI. RESULTS AND DISCUSSION In this section we employ a simple model system, i.e., a dissipative two-site system, to validate the Chebyshev QME and to compare the results to the aforementioned established methods. To this end, we consider the system Hamiltonian of form discussed before, eq 1, with site energies E1 and E2 and coupling elements V12 = V21 = V. Each site is coupled independently to its own thermal bath, where the coupling is governed by Drude−Lorentz spectral densities defined as ω Ji (ω) = 2πλiγi 2 ω + γi 2 (18) with the real parameters γi and λi. Specifically, γi affects the width of the function and represents the inverse correlation time, whereas the reorganization energy λi influences its amplitude. For the sake of simplicity, however, we chose the same spectral density for both sites. This system has been used in several benchmark studies in literature41−43 in which the 3274
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A of equal site energies. Again, the CDTL scheme almost reaches the same equilibrium value as the converged results. Small deviations in the thermal state can be explained by the secondorder nature of the present approach.45 In addition, the Markovian limit of the present scheme is shown in Figures 1 and 2. The procedure how this limit has been obtained, was detailed in section IV. For the case of equal site energies, the Markovian limit actually shows larger oscillations than the non-Markovian results and is less accurate at the same time. This behavior is different for the case of unequal site energies. Not many differences are visible for the smallest reorganization energy. For the reorganization energy of λ = 20 cm−1 the agreement with the HEOM result is better for the Markovian than for the non-Markovian dynamics. For even larger reorganization energies, the reminiscences of oscillations, present in the non-Markovian case, vanish completely in the Markovian limit. This behavior of the Markovian case agrees well with the Markovian Redfield dynamics as reported earlier.41,46 Moreover, shown in Figures 1 and 2 are the results for the classical path-based ensemble-average wave packet dynamics outlined in the previous section. It is interesting to see how different this approximation behaves compared to a perturbative treatment. For the case of equal site energies in Figure 1, the wave packet dynamics adequately agrees with the HEOM results even for reorganization energies up to λ = 100 cm−1. Only for the largest reorganization energy quantitative differences are visible, though qualitatively the behavior is similar again. So for equal site energies, the ensemble-average wave packet dynamics approach outperforms the present Chebyshev scheme, which is based on a second-order treatment. In the wave packet dynamics approach the dephasing is actually included properly in the theory by the ensemble average. This theory is missing a relaxation term, however. A fact that is invisible for equal site energies. Moving on to a bias between the site energies this deficiency of the wave packet dynamics approach becomes obvious because the long-term limit is not properly reproduced, a fact that is wellknown.43 Nevertheless, it is clearly visible that the dephasing part of the dynamics is still quite accurately reproduced, although the dynamics lead to the wrong thermal state. Finally, in Figure 3 we deal with a nonsmooth spectral density as depicted in the inset. The spectral density linearly increases with the frequency and then suddenly drops to zero. Such a spectral density is practically not reproducible by a sum of Lorentzian functions; i.e., it cannot be handled properly either by the EDTL or by standard HEOM approaches. Note that the spectral density is extended to negative frequencies by Ji(−ω) = −Ji(ω) before being inserted into eq 5. The thermal equilibria are the same as in Figures 1 and 2, respectively, because these do not depend on the spectral density when higher-order corrections are neglected.45 This example clearly shows the advantage of the present scheme in being able to use any form of spectral density. Similar examples can be tested for low temperatures though one runs the risk of leaving the validity range of the second-order perturbation theory.
Figure 3. Population dynamics of the first site determined by the CDTL scheme, for both configurations, i.e., ΔE = 0 and ΔE = 100 cm−1. The employed spectral density is shown in the inset. It consists of a linearly increasing function which drops to zero at half the maximum frequency. The maximum of the spectral density Jmax is chosen such that it should lie within the validity range of perturbation theory.
reservoirs.25 As in earlier approaches,14,23,24 including the EDTL one,22 the present scheme does not make any additional assumptions and approximations in the time dependence of the system Hamiltonian. So the approach is able to treat systems with rapidly fluctuating and/or strong external or internal perturbations. In the Fermionic case this has been employed, for example, in studying the coherent destruction of tunneling in molecular wire systems.47 In these investigations we used the EDTL scheme which is only efficiently usable when the spectral density can be decomposed in terms of a few Lorentzian functions and if the temperature is not too low because otherwise too many terms are needed in the decomposition of the Fermi function. The CDTL approach lifts the restriction of having to use simple forms of spectral densities and not too low temperatures. The drawback is that the simulation time is limited by the number of terms used in the Chebyshev expansion. As detailed in eq 15, the present scheme depends quasi-linearly on the simulation time; i.e., when the damping is small resulting in a longer simulation time needed until the steady state is reached, the scheme becomes less efficient with regard to computational time. At the same time, the proposed scheme depends quasi-linearly on the spectral range Ω of the underlying system as well. So it becomes more efficient for systems with a narrow spectral range. Moreover, a slight dependence on the reorganization energy of the environmental coupling was observed in the numerical studies which needs to be analyzed in more detail in the future. In recent years complex spectral densities have been extracted from experiments48,49 or from a combination of classical molecular dynamics simulations and quantum chemical calculations.20,21,50−53 These spectral densities show a very complex structure which can only approximately be reproduced by a sum of Lorentzian functions.8 To be able to treat systems with these spectral densities, one probably first has to extend the present scheme to higher orders or to the HEOM scheme, but we have clearly shown that such an approach based on a Chebyshev decomposition is feasible and a practical alternative to the exponential decomposition schemes.
VII. CONCLUDING REMARKS In this contribution we have introduced a new quantum master equation based on a second-order treatment in the environmental coupling employing a Chebyshev decomposition of the bath correlation function. The CDTL scheme is an extension of a recently reported approach for coupling to Fermionic 3275
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A
■
(18) Xie, H.; Kwok, Y.; Zhang, Y.; Jiang, F.; Zheng, X.; Yan, Y.; Chen, G. Time-dependent Quantum Transport Theory and Its Applications to Graphene Nanoribbons. Phys. Status Solidi B 2013, 250, 2481. (19) Ritschel, G.; Eisfeld, A. Analytic Representation of Bath Correlation Functions for Ohmic and Superohmic Spectral Densities Using Simple Poles. J. Chem. Phys. 2014, 141, 094101. (20) Olbrich, C.; Kleinekathöfer, U. Time-dependent Atomistic View on the Electronic Relaxation in Light-harvesting System II. J. Phys. Chem. B 2010, 114, 12427−12437. (21) Jurinovich, S.; Viani, L.; Curutchet, C.; Mennucci, B. Limits and Potentials of Quantum Chemical Methods in Modelling Photosynthetic Antennae. Phys. Chem. Chem. Phys. 2015, 17, 30783−30792. (22) Kleinekathöfer, U. Non-Markovian Theories Based on the Decomposition of the Spectral Density. J. Chem. Phys. 2004, 121, 2505. (23) Welack, S.; Schreiber, M.; Kleinekathöfer, U. The Influence of Ultra-fast Laser Pulses on Electron Transfer in Molecular Wires Studied by a Non-Markovian Density Matrix Approach. J. Chem. Phys. 2006, 124, 044712. (24) Tanimura, Y. Stochastic Liouville, Langevin, Fokker-Planck, and Master Equation Approaches to Quantum Dissipative Systems. J. Phys. Soc. Jpn. 2006, 75, 082001. (25) Popescu, B.; Rahman, H.; Kleinekathöfer, U. Using the Chebyshev Expansion in Quantum Transport Calculations. J. Chem. Phys. 2015, 142, 154103. (26) Kosloff, R. Time-dependent Quantum-mechanical Methods for Molecular Dynamics. J. Phys. Chem. 1988, 92, 2087−2100. (27) Kosloff, R. Propagation Methods for Quantum Molecular Dynamics. Annu. Rev. Phys. Chem. 1994, 45, 145. (28) Chen, R.; Guo, H. The Chebyshev Propagator for Quantum Systems. Comput. Phys. Commun. 1999, 119, 19. (29) Vijay, A.; Metiu, H. A Polynomial Expansion of the Quantum Propagator, the Green’s Function, and the Spectral Density Operator. J. Chem. Phys. 2002, 116, 60−68. (30) Ndong, M.; Tal Ezer, H.; Kosloff, R.; Koch, C. P. A Chebychev Propagator for Inhomogeneous Schrödinger Equations. J. Chem. Phys. 2009, 130, 124108. (31) Ndong, M.; Tal Ezer, H.; Kosloff, R.; Koch, C. P. A Chebychev Propagator with Iterative Time Ordering for Explicitly Timedependent Hamiltonians. J. Chem. Phys. 2010, 132, 064105. (32) Guo, H.; Chen, R. Short-time Chebyshev Propagator for the Liouville-von Neumann Equation. J. Chem. Phys. 1999, 110, 6626. (33) Huisinga, W.; Pesce, L.; Kosloff, R.; Saalfrank, P. Faber and Newton Polynomial Integrators for Open-system Density Matrix Propagation. J. Chem. Phys. 1999, 110, 5538. (34) Woods, M. P.; Groux, R.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Mappings of Open Quantum Systems onto Chain Representations and Markovian Embeddings. J. Math. Phys. 2014, 55, 032101. (35) Arfken, G. Mathematical Methods for Physicists, 3rd ed.; Academic Press: San Diego, 1985. (36) de Doncker-Kapenga, E.; Kahaner, D. K.; Piessens, R.; Ü berhuber, C. W. Quadpack: A Subroutine Package for Automatic Integration; Springer: Berlin, 1983. (37) Coifman, R.; Rokhlin, V.; Wandzura, S. The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription. Antennas and Propagation Magazine, IEEE 1993, 35, 7−12. (38) Landa, Y.; Tanushev, N.; Tsai, R. Discovery of Point Sources in the Helmholtz Equation Posed in Unknown Domains with Obstacles. Comm. in Math. Sci. 2011, 9, 903−928. (39) Ishizaki, A.; Tanimura, Y. Quantum Dynamics of Systems Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierachy Equations Approach. J. Phys. Soc. Jpn. 2005, 74, 3131−3134. (40) Amin, A. F.; Li, G.-Q.; Phillips, A. H.; Kleinekathöfer, U. Coherent Control of the Spin Current Through a Quantum Dot. Eur. Phys. J. B 2009, 68, 103−109. (41) Ishizaki, A.; Calhoun, T. R.; Schlau Cohen, G. S.; Fleming, G. R. Quantum Coherence and Its Interplay with Protein Environments in Photosynthetic Electronic Energy Transfer. Phys. Chem. Chem. Phys. 2010, 12, 7319−7337.
AUTHOR INFORMATION
Corresponding Author
*U. Kleinekathö f er. E-mail:
[email protected]. Present Address ‡
Max-Planck-Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187, Dresden, Germany. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Financial support by the Deutsche Forschungsgemeinschaft (DFG) through project KL1299/14-1 is gratefully acknowledged.
■
REFERENCES
(1) Weiss, U. Quantum Dissipative Systems, 2nd ed.; World Scientific: Singapore, 1999. (2) Breuer, H. P.; Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press: Oxford, U.K., 2002. (3) May, V.; Kühn, O. Charge and Energy Transfer in Molecular Systems, 3rd ed.; Wiley-VCH: Berlin, 2011. (4) Tanimura, Y.; Kubo, R. Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath. J. Phys. Soc. Jpn. 1989, 58, 101. (5) Jin, J. S.; Zheng, X.; Yan, Y. J. Exact Dynamics of Dissipative Electronic Systems and Quantum Transport: Hierarchical Equations of Motion Approach. J. Chem. Phys. 2008, 128, 234703. (6) Ishizaki, A.; Fleming, G. R. Theoretical Examination of Quantum Coherence in a Photosynthetic System at Physiological Temperature. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 17255−17260. (7) Strü m pfer, J.; Schulten, K. Open Quantum Dynamics Calculations with the Hierarchy Equations of Motion on Parallel Computers. J. Chem. Theory Comput. 2012, 8, 2808−2816. (8) Kreisbeck, C.; Kramer, T.; Aspuru-Guzik, A. Scalable HighPerformance Algorithm for the Simulation of Exciton Dynamics. Application to the Light-Harvesting Complex II in the Presence of Resonant Vibrational Modes. J. Chem. Theory Comput. 2014, 10, 4045−4054. (9) Schrö der, M.; Schreiber, M.; Kleinekathö fer, U. Reduced Dynamics of Coupled Harmonic and Anharmonic Oscillators Using Higher-order Perturbation Theory. J. Chem. Phys. 2007, 126, 114102. (10) Strümpfer, J.; Schulten, K. Excited State Dynamics in Photosynthetic Reaction Center and Light Harvesting Complex 1. J. Chem. Phys. 2012, 137, 065101−8. (11) Tian, H.; Chen, G. An Efficient Solution of Liouville-von Neumann Equation That Is Applicable to Zero and Finite Temperatures. J. Chem. Phys. 2012, 137, 204114. (12) Tian, H.; Chen, G. Application of Hierarchical Equations of Motion (HEOM) to Time Dependent Quantum Transport at Zero and Finite Temperatures. Eur. Phys. J. B 2013, 86, 411. (13) Popescu, B.; Kleinekathöfer, U. Treatment of Time-dependent Effects in Molecular Junctions. Phys. Status Solidi B 2013, 250, 2288− 2297. Special issue ’Quantum Transport at the Molecular Scale’ (14) Meier, C.; Tannor, D. J. Non-Markovian Evolution of the Density Operator in the Presence of Strong Laser Fields. J. Chem. Phys. 1999, 111, 3365. (15) Hu, J.; Xu, R. X.; Yan, Y. J. Pade Spectrum Decomposition of Fermi Function and Bose Function. J. Chem. Phys. 2010, 133, 101106. (16) Croy, A.; Saalmann, U. Partial Fraction Decomposition of the Fermi Function. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 159904. (17) Liu, H.; Zhu, L.; Bai, S.; Shi, Q. Reduced Quantum Dynamics with Arbitrary Bath Spectral Densities: Hierarchical Equations of Motion Based on Several Different Bath Decomposition Schemes. J. Chem. Phys. 2014, 140, 134106. 3276
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
Article
The Journal of Physical Chemistry A (42) Ishizaki, A.; Fleming, G. R. On the Interpretation of Quantum Coherent Beats Observed in Two-Dimensional Electronic Spectra of Photosynthetic Light Harvesting Complexes. J. Phys. Chem. B 2011, 115, 6227−6233. (43) Aghtar, M.; Liebers, J.; Strü m pfer, J.; Schulten, K.; Kleinekathöfer, U. Juxtaposing Density Matrix and Classical Pathbased Wave Packet Dynamics. J. Chem. Phys. 2012, 136, 214101. (44) van der Vegte, C. P.; Dijkstra, A. G.; Knoester, J.; Jansen, T. L. C. Calculating Two-Dimensional Spectra with the Mixed QuantumClassical Ehrenfest Method. J. Phys. Chem. A 2013, 117, 5970−5980. (45) Geva, E.; Rosenman, E.; Tannor, D. J. On the Second-order Corrections to the Quantum Canonical Equilibrium Density Matrix. J. Chem. Phys. 2000, 113, 1380. (46) Ishizaki, A.; Fleming, G. R. Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach. J. Chem. Phys. 2009, 130, 234111. (47) Li, G.-Q.; Schreiber, M.; Kleinekathöfer, U. Coherent Laser Control of the Current Through Molecular Junctions. EPL 2007, 79, 27006. (48) Adolphs, J.; Renger, T. How Proteins Trigger Excitation Energy Transfer in the Fmo Complex of Green Sulfur Bacteria. Biophys. J. 2006, 91, 2778−2797. (49) Kell, A.; Feng, X.; Reppert, M.; Jankowiak, R. On the Shape of the Phonon Spectral Density in Photosynthetic Complexes. J. Phys. Chem. B 2013, 117, 7317−7323. (50) Valleau, S.; Eisfeld, A.; Aspuru Guzik, A. On the Alternatives for Bath Correlators and Spectral Densities from Mixed QuantumClassical Simulations. J. Chem. Phys. 2012, 137, 224103−13. (51) Aghtar, M.; Strü m pfer, J.; Olbrich, C.; Schulten, K.; Kleinekathöfer, U. Different Types of Vibrations Interacting with Electronic Excitations in Phycoerythrin 545 and Fenna-MatthewsOlson Antenna Systems. J. Phys. Chem. Lett. 2014, 5, 3131−3137. (52) Viani, L.; Corbella, M.; Curutchet, C.; O’Reilly, E. J.; Olaya Castro, A.; Mennucci, B. Molecular Basis of the Exciton-phonon Interactions in the PE545 Light-harvesting Complex. Phys. Chem. Chem. Phys. 2014, 16, 16302−16311. (53) Jurinovich, S.; Curutchet, C.; Mennucci, B. The FennaMatthews-Olson Protein Revisited: A Fully Polarizable (TD)DFT/ MM Description. ChemPhysChem 2014, 15, 3194−3204.
3277
DOI: 10.1021/acs.jpca.5b12237 J. Phys. Chem. A 2016, 120, 3270−3277
