Iterative Calculation of Energy Eigenstates Employing the Multilayer

May 15, 2014 - A large number of degrees of freedom, up to 10 000 bath modes, can be treated, and thus the appropriate phase transition at a critical ...
0 downloads 0 Views 618KB Size
Article pubs.acs.org/JPCA

Iterative Calculation of Energy Eigenstates Employing the Multilayer Multiconfiguration Time-Dependent Hartree Theory Haobin Wang* Beijing Computational Science Research Center, No. 3 He-Qing Road, Hai-Dian District, Beijing 100084, P.R. China Department of Chemistry and Biochemistry, MSC 3C, New Mexico State University, Las Cruces, New Mexico 88003, United States ABSTRACT: The improved relaxation approach of Meyer and co-workers is extended by employing the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) theory. Within this approach the trial wave function is expressed by a recursive, multilayered expansion that is the same as in a normal ML-MCTDH simulation. Calculation of an energy eigenstate then proceeds by iteratively diagonalizing the Boltzmann operator in the top layer representation using a Lanczos/Arnoldi method while relaxing the single particle functions of all layers using the ML-MCTDH imaginary time propagation. Numerical iteration is stopped until the relative energy and/or density of the target state meet the convergence criteria. The application of this multilayer improved relaxation method is illustrated by computing the energy splitting for the lowest pair of eigenstates as well as the spin observables and entanglement entropy of the spin-boson model. A large number of degrees of freedom, up to 10 000 bath modes, can be treated, and thus the appropriate phase transition at a critical system-bath coupling strength may be predicted via appropriate extrapolations. systems.31−39 Within this approach the wave function is expressed by a recursive, layered expansion,

I. INTRODUCTION The description of many-body quantum dynamical effects for reactions in complex systems is a challenging task in theoretical physical chemistry. According to the nature of the methods the current development can be roughly divided into two classes, approximate and numerically exact. When treating many physical models, quantum perturbation theories are often used to describe the reduced dynamics of a subsystem.1−6 In chemistry where numerical simulation is often employed, the popular approximations are mixed quantum-classical and semiclassical approaches.7−13 Such methods are applicable with analytic potential functions or even ab initio potential energy surfaces generated on-the-fly and are thus becoming increasingly more widespread. Another important direction is the development of numerically exact approaches, e.g., numerical path integral approaches14−18 using Feynman−Vernon influence functional,19 the numerical renormalization group (NRG) theory,20−22 and the numerically exact self-consistent hybrid method.23−25 Complementary to the approximate approaches, these methods aim at providing accurate results for physically motivated models whose parameters may be obtained from electronic structure calculations and/or experiments. The accurate simulation on these model systems may be regarded as “numerical experiments” with reliable results, which is often the first step in developing a general theory for describing the underlying processes. Along this line we have focused on developing the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) theory,26 which was motivated by the original MCTDH method27−30 and extends its applicability to significantly larger © 2014 American Chemical Society

p

|Ψ(t )⟩ =

∑ ∑ ··· ∑ A j j ...j (t ) ∏ |φ j(κ)(t )⟩ 12

j1

|φ j(κ)(t )⟩ = κ

j2

p

κ=1

jp

∑ ∑ ··· ∑ Biκi,j...i κ

12

i1

i2

iQ (κ)

(t ) Q (κ )

κ

(1.1a) Q (κ )

∏ |vi(κ ,q)(t )⟩ q

q=1

(1.1b) M(κ , q)

|vi(qκ , q)(t )⟩ =

∑ ∑ ··· ∑ α1

α2

αM(κ , q)

Cακ1,αq2,...iq αM(κ ,q)(t )



|ξα(γκ , q , γ )(t )⟩

γ=1

(1.1c)

... κ,jκ κ where Aj1j2...jp(t), Bκ,j i1i2...iQ(κ)(t), Cα1α2...αM(κ,q)(t), and so on are the expansion coefficients for the first, second, third, ..., layers, (κ,q) (κ,q,γ) respectively; |φ(κ) (t)⟩, ..., are the “single jκ (t)⟩, |viq (t)⟩, |ξαγ particle” functions (SPFs) for the first, second, third, ..., layers. The multilayer hierarchy is terminated at a particular level by expanding the SPFs in the deepest layer in terms of timeindependent configurations, each of which may contain several Cartesian degrees of freedom. The variational parameters

Special Issue: International Conference on Theoretical and High Performance Computational Chemistry Symposium Received: April 4, 2014 Revised: May 10, 2014 Published: May 15, 2014 9253

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

within the ML-MCTDH theoretical framework are dynamically optimized through the use of Dirac−Frenkel variational principle,40 ⟨δΨ(t)|i(∂/∂t) − Ĥ |Ψ(t)⟩ = 0, which results in a set of coupled, nonlinear differential equations for the expansion coefficients of all layers.26,37,41 The theory has also been generalized to treat quantum identical particles explicitly by introducing the multilayer variation in the abstract Fock space using the occupation number representation, i.e., the MLMCTDH theory in the second quantized representation.41 This unified ML-MCTDH theory has been developed over the past few years and has seen promising applications.42−46 The introduction of the recursive, dynamically optimized layering scheme in the ML-MCTDH wave function provides more flexibility in the variational functional, which results in a tremendous gain in one’s ability to study large many-body quantum systems. This is demonstrated by many applications on simulating quantum dynamics of ultrafast electron transfer reactions in condensed phases.32−34,36−39,47−54 The MLMCTDH work of Manthe has introduced an even more adaptive formulation based on a layered correlation discrete variable representation.55,56 The applicability of the ML-MCTDH theory is not only limited to simulating quantum dynamics. Similar to the previous formulation of MCTDH for calculating eigenstates of the Hamiltonian30,57,58 or the thermal flux operator,58,59 the ML-MCTDH theory can be applied in the same way. This extension may be useful in several aspects. For example, there are situations where a correlated initial state is required, e.g., simulating time-resolved spectroscopy in which the material system is initially at thermal equilibrium, calculating the reactive flux correlation function, or studying electron transport with an initial grand canonical ensemble at zero bias voltage. At zero (or a very low) temperature, a pure relaxation approach, based entirely on imaginary time propagation, may be prohibitively expensive. Thus, a more efficient method is desirable. Another application would be the calculation of a few lowest energy eigenstates under critical conditions so as to predict important thermodynamic properties of a large system, e.g., the quantum phase transition of the spin-boson model to be illustrated later. In this paper we discuss the extension of the ML-MCTDH theory to facilitate improved relaxation simulation proposed by Meyer and co-workers.57 Section II outlines both the formal theory and practical implementation within the ML-MCTDH framework. Section III illustrates the application of this approach to the calculation of energy splitting and spin observables for the spin-boson model in various parameter regimes, and how the result can be used to predict the quantum phase transition in this model. Section IV concludes.

δ{⟨Ψ|Ĥ |Ψ⟩ − E(⟨Ψ|Ψ⟩ − 1) p



∑ ∑ ϵ(jκ,l)(⟨φj(κ)|φl(κ)⟩ − δjl) κ=1 j,l p



Q (κ )

∑ ∑ ∑ ϵ̃(jκ,l,q)(⟨v(j κ ,q)|vl(κ ,q)⟩ − δjl) κ=1 q=1 p



j,l

Q (κ ) M(κ , q)

∑ ∑ ∑ ∑ ϵ̃∼ j , l κ=1 q=1

γ=1

(κ , q , γ )

(⟨ξj(κ , q , γ )|ξl(κ , q , γ )⟩ − δjl) − ···}

j,l

(2.1)

=0

where E is the Lagrange multiplier to keep the overall wave (κ,q) ∼ (k,q,γ) , etc. are Lagrange function normalized and ϵ(κ) j,l , ϵ̃j,l , ϵ̃j,l multipliers to keep the SPFs of the first, second, third, etc. layers orthonormal. To carry out the variation in (2.1), we assume that the expansion coefficients of all the layers, e.g., κ κ,q,iq Aj1j2...jp,, Bκ,j i1i2...iQ(κ) Cα1α2...αM(κ,q), etc. in (1.1), are complex, so that their complex conjugates (bra) are varied independently (if they are real, the final result is the same.) Variation in the top layer gives the standard eigen-equation HJLAL ≡

∑ ⟨ΦJ |Ĥ |ΦL⟩AL = EAJ L

(2.2)

where we have denoted a configuration |ΦJ⟩ as the Hartree product |ΦJ⟩ ≡ ∏pκ=1|φ(κ) jκ and used the collective index, Aj1j2...jp ≡ AJ. Solving this eigen-equation thus gives the eigen-energy E and the eigenvector Ψ = ∑JAJ|ΦJ⟩ in terms of the top layer expansion coefficients AJ and the SPFs of all the layers (because Ψ = ∑JAJ|ΦJ⟩ depends on them). Variation (2.1) with respect to the SPFs of all the layers, after eliminating all the Lagrange multipliers, have the same form (κ ) [1 − P ̂ ]⟨Ĥ ⟩(κ)| φ̲ (κ)⟩ = 0

(2.3a)

(κ , q) [1 − PL̂ 2 ]⟨/̂ ⟩(L2κ , q)| v(̲ κ , q)⟩ = 0

(2.3b)

(κ , q , γ )

[1 − PL̂ 3

]⟨/̂ ⟩(L3κ , q , γ )| ξ̲ (κ , q , γ )⟩ = 0

(2.3c)

where ⟨Ĥ ⟩(κ), ⟨/̂ ⟩(L2κ , q), ⟨/̂ ⟩(L3κ , q , γ ), and so on are mean-field operators for the first, second, third, and so on layers, ̂(κ,q,γ) respectively; P̂(κ), P̂ (κ,q) L2 , PL3 , and so on are projection operators in the single particle space using the SPFs of the T (κ,q) (κ) corresponding layers; and |φ(k)⟩ = {|φ(κ) ⟩= 1 ⟩, |φ2 ⟩, ...} , |v T (κ,q,γ) T (κ,q) (κ,q) (κ,q,γ) (κ,q,γ) {|v1 ⟩, |v2 ⟩, ...} , |ξ ⟩ = {|ξ1 ⟩, |ξ2 ⟩, ...} , etc. denote the symbolic column vector of the SPFs for different layers. The definition of these intermediate quantities uses the SPFs or the single hole functions as detailed in the previous MLMCTDH publications.26,37,41 One may propose to solve the above set of equations, eqs 2.2−2.3, iteratively so that eventually the eigen-energy E and the eigenvector Ψ = ∑JAJ|ΦJ⟩ can be obtained. However, just as in the early days of MCSCF theory,60 the resulting equations are highly nonlinear to iterate and very difficult to converge. To get around this difficulty, the MCTDH “improved relaxation” method57 uses repeated, short imaginary time propagation to relax the SPFs, which is then combined with a Krylov subspace method to solve the eigen-equation (2.2). This is generalized in a straightforward way employing the ML-MCTDH method. B. Practical Implementation. Similar to the single layer MCTDH case as described in ref 57, the stationary solution in

II. MULTILAYER IMPROVED RELAXATION A. Formal Theory. To compute an energy eigenstate with the multilayer wave function ansatz (except that the coefficients are no longer time-dependent), eq 1.1, one applies the timeindependent variation analogous to that used in ref 57 9254

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

One may also use Davidson’s method as in the previous work by Meyer and co-workers57 with a preconditioner, e.g., Jacobi (diagonal), Gauss−Seidel, or symmetric successive over relaxation (SSOR)62 to accelerate the Krylov iteration. After the top layer Hamiltonian matrix is diagonalized, the Acoefficients of the target eigenstate are fixed while the SPFs of all layers are relaxed over a certain imaginary time interval by solving eq 2.5. This completes one step of iteration. In the next step, the Hamiltonian matrix at the top layer is rebuilt using the relaxed SPFs and diagonalized again, followed by the imaginary time relaxation of the SPFs. This process is iterated until convergence, usually defined by the relative energy or density between two successive iterations, is achieved for the target state of interest. For a large system with many degrees of freedom, it is practically impossible to solve for all the eigenstates. However, there is often physical significance to analyze states with the lowest energies, e.g., the ground state and some low-lying excited states. As discussed in the Introduction, such analysis may be useful for both thermodynamic study and quantum dynamics simulation at low temperatures where correlated initial states are required. This makes the multilayer improved relaxation approach a valuable tool. As an example, we will apply our method to study the spin-boson model.

(2.3) remains valid under any linear transformation. Thus, one may insert the pseudoinverse of the reduced density matrices (defined recursively in terms of the hole functions26,37,41) ρ̂(κ), (κ,q,γ) ϱ̂(κ,q) L2 , ϱ̂ L2 , ..., to the equation of the corresponding layer (κ ) [1 − P ̂ ][ρ(̂ κ) ]−1 ⟨Ĥ ⟩(κ)⟨ φ̲ (κ)⟩ = 0

(2.4a)

(κ , q)

̂ ][ϱ̂(κ , q)]−1 ⟨/̂ ⟩(L2κ , q)⟨ v(̲ κ , q)⟩ = 0 [1 − PL2 L2

(2.4b)

(κ , q , γ ) ̂ [1 − PL3 ][ϱ̂(L3κ , q , γ )]−1 ⟨/̂ ⟩(L3κ , q , γ )| ξ̲ (κ , q , γ )⟩ = 0

(2.4c)

and realize that these are related to the imaginary time derivative of the SPFs, i.e. (κ ) φ̲ ̇k (τ )⟩L2coefficients = −[1 − P ̂ ][ρ(̂ κ) ]−1 ⟨Ĥ ⟩(κ)| φ̲ (κ)⟩

(2.5a)

| v̲ (̇ κ , q)(τ )⟩L3coefficients (κ , q) ̂ ][ϱ̂(κ , q)]−1 ⟨/̂ ⟩(L2κ , q)| v(̲ κ , q)⟩ = −[1 − PL2 L2

(2.5b)

(κ , q , γ )

| ξ̲ ̇

(τ )⟩L4coefficients (κ , q , γ )

̂ = −[1 − PL3

][ϱ̂(L3κ , q , γ )]−1 ⟨/̂ ⟩(L3κ , q , γ )| ξ̲ (κ , q , γ )⟩

(2.5c)

III. ILLUSTRATIVE EXAMPLE: THE SPIN-BOSON MODEL A. Overview. The spin-boson Hamiltonian3,63 is widely used in condensed phase physics and chemistry to model a variety of different processes such as electron transfer,64 hydrogen tunneling,65 macroscopic quantum coherence,66 and many others.63 The model comprises two electronic states linearly coupled to a bath of harmonic oscillators. In massweighted coordinates the Hamiltonian reads

For clarity we refer to the top layer as the level one (L1) SP space, the second layer as the level 2 (L2) SP space, and so on. The SPF of a particular layer depends the expansion coefficients of the layer below, in the form of eq (1.1). In the notation used in eq (2.5) it is implicitly understood that the imaginary time derivatives on the left-hand side are only performed with respect to the expansion coefficients of a particular layer (denoted by the respective subscript). For example, the time derivative in eq 2.5a is on the L2 expansion κ coefficient Bκ,j i1,i2...iQ(κ)(t), etc. Finally, the imaginary time is defined as τ = −it. Thus, the approximate solution to eq (2.3) is to use the MLMCTDH method to propagate the SPFs in imaginary time, applying eq 2.5, until the time derivatives are sufficiently small. This is then coupled to the solution of the top layer eigenequation (2.2) in an iterative way. Overall, the multilayer improved relaxation method employs a procedure similar to that discussed in ref 57. An initial guess is usually generated by diagonalizing a zeroth-order Hamiltonian Ĥ 0 and choosing an appropriate target (ground or excited) state. Besides the SPFs that are used to define this initial state, other SPFs are chosen randomly and orthogonalized. This constructs a multilayer wave function in the form of eq 1.1, often with only a few nonzero top layer A-coefficients. The wave function is relaxed over a short imaginary time to generate a correlated initial state, where most A-coefficients are no longer zero. Then one proceeds to the iterations of multilayer improved relaxation as follows. In each iteration step, the Hamiltonian matrix of the top layer is first constructed on the basis of the current SPFs. Then one proceeds to diagonalize this Hermitian matrix as shown by eq 2.2. Due to the fact that the matrix is often big in size and in the mean time one is often only interested in a few eigenstates with extreme eigenvalues, an iterative solver such as a Krylov subspace method is usually the choice. In this work an Arnoldi/ Lanczos method on the Boltzmann operator is employed, in spirit similar to the work of Manthe58,61 but only applying the Krylov iteration to the linear part of the top layer A-coefficients.

H = ϵσz + Δσx +

1 2

∑ (pi 2

+ ωi 2qi 2) + σz ∑ ciqi

i

i

(3.1)

where σx and σz are the Pauli matrices σx = |ϕ1⟩⟨ϕ2| + |ϕ2⟩⟨ϕ1|

(3.2a)

σz = |ϕ1⟩⟨ϕ1| − |ϕ2⟩⟨ϕ2|

(3.2b)

The properties of the bath that influence the dynamics of the two-state subsystem are completely specified by the spectral density function3,63 JB (ω) =

π 2

∑ j

cj 2 ωj

δ(ω − ωj) (3.3)

Despite its simple form, the spin-boson problem appears to be challenging with respect to both analytical and numerical treatments. This is especially the case when the properties in the zero temperature limit of the phonon bath are studied, which is characterized by strong quantum effects. Apart from the insight offered from time-dependent studies,36,37 it is desirable to analyze the infinite time thermodynamic limit by a very accurate method. This is particularly useful when one investigates the effect of destructive quantum interference, that is, phonon-induced localization at zero temperature.3,20,63,67 According to previous studies, the quantum tunneling between the two states in the spin-boson model disappears when the system-bath coupling goes beyond a critical strength. 9255

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

according to the number of the bath modes, e.g., two layers for a bath of 100−300 modes, three layers for a bath of 400−1000 modes, and four layers for a bath up to 10000 modes. After the number of layers is determined, the bath modes are lined up according to the order of ascending/descending frequencies and divided into different groups/subgroups by joining the nearest neighbors. Eventually all the degrees of freedom are arranged under a multilayer tree. The structure of the spinboson Hamiltonian makes it easy to do so. This uniform, balanced tree structure is also due to the fact that the bath discretization scheme described above ensures that the coupling strength of each bath mode is roughly the same. In situations where many strongly coupled intramolecular modes exist with different coupling strength,50 a more delicate multilayer tree set up would be required, usually involving some empirical adjustment. Let us consider an example where a three-layer ML-MCTDH scheme is used to treat a bath of 1000 modes. In this case a tree structure is set up as follows: (1) In the top layer (L1 layer) there are five (four bath plus one electronic) composite L1 single particle groups. (2) In the second layer (L2 layer) there are four L2 single particle subgroups within each bath L1 single particle group of the previous layer, making a total of 4 × 4 = 16 single particle groups in the L2 layer. (3) In the third layer (L3 layer) there are again four L3 single particle subgroups chosen within each L2 single particle group, resulting in a total of 4 × 4 × 4 = 64 single particle groups in the L3 layer. This completes the overall multilayer tree structure for the calculation. The deepest layer then employs static basis functions. For each mode the primitive basis functions are the eigenfunctions of this harmonic oscillator. Naturally, if a mode is anharmonic, eigenfunctions of the corresponding anharmonic mode is a sound choice.34 The number of the primitive basis functions is chosen according to a maximum cutoff energy (or, given a temperature, a smallest cutoff probability of occupation). In the present work this number is between three and a few hundred depending on the frequency of the mode. Because ML-MCTDH is not as restricted by the number of basis functions as some other methods, we often use more than enough primitive basis functions. At this point, the task is to distribute 1000 modes evenly among the 64 L3 single particle groups. The program uses an iterative procedure to achieve this. First, an initial partition of the 1000 modes is generated where 12−20 modes are allocated to each L3 single particle group on the basis of any reasonable guess (“reasonable” in a sense that more modes are combined in one group at the higher frequency end of the spectral density because they have fewer basis functions). The number of the primitive basis functions within each L3 single particle group is the size of the full configurations within this group and thus a huge number. These primitive basis functions are adiabatically contracted26,37 according to the criterion that the total zerothorder energy within each L3 group is below a maximum energy (in our program this is a factor of 2−3 of the cutoff energy for a single mode, which is a convergence parameter). After the contraction, the numbers of (contracted) basis functions between different L3 single particle groups are often quite different. To balance the number of basis sets, an improved partition scheme is generated where modes from a group with bigger contracted basis size are moved to a smaller basis size group. Adiabatic basis set contraction is carried out again to these newly generated L3 single particle groups. The procedure

This quantum phase transition could not be obtained from a perturbative approach but has been predicted by various analytical theories,3,63,68,69 formal mathematical arguments,67 the flow-equation approach70 as well as NRG calculations.20,21 In the present work we use a procedure similar to that described in a previous study employing the density matrix renormalization group (DMRG) approach,71 where the energy splitting for the lowest pair of eigenstates are examined versus the number of phonon modes. Because our approach also generates the wave function, any other physical quantities may also be evaluated. For example, the von Neumann entanglement entropy and the longitudinal spin magnetization ⟨σz⟩ discussed in the previous work.72 B. Details of the Calculation. As discussed above, the bath and its coupling to the two-level system is fully characterized by the spectral density JB(ω) in eq 3.3. In this work we consider a spectral density of Ohmic form with an exponential cutoff π J(ω) = αωe−ω / ωc (3.4) 2 Here, α is the dimensionless Kondo parameter that characterizes the system−bath coupling strength and ωc is the characteristic frequency of the bath. The relation of the spinboson model to the Kondo problem has been discussed previously.3,63 Other more complex forms of the spectral density or anharmonic bath have been studied previously at finite temperatures.34,38,39,73 The thermodynamics of these models is the subject of future studies. The continuous bath spectral density of eq 3.4 can be discretized to the form of eq 3.3 via the relation23,24,26 cj 2 =

2 J(ωj) ωj π ρ(ωj)

(3.5)

where ρ(ω) is a density of frequencies satisfying

∫0

ωj

dω ρ(ω) = j

j = 1, ..., Nb

(3.6a)

In this work, ρ(ω) is chosen as ρ(ω) ∝

N + 1 −ω / ωc J(ω) = b e ω ωc

(3.6b)

where Nb is the number of the discrete bath modes in the simulation. We note that if there is no exponential damping term in eq 3.3, i.e., a strict linear (Ohmic) expression, and then ρ(ω) becomes a constant. The discretization gives equal spacing in frequencies between adjacent harmonic modes as in the DMRG work.71 For the example considered in this paper we wish to examine the dependence of a physical quantity on the number of the bath modes Nb. This way a thermodynamic limit may be extracted in the limit of an infinite bath. Each fixed Nb defines a (large) finite system, for which multilayer improved relaxation calculations are carried out with all the variational parameters systematically increased to achieve convergence. In this study, we have used the criterion that the relative error for the energy of the target state be smaller than a preset tolerance (10−8− 10−6) within each relaxation calculation. When the numbers of basis functions, configurations, etc. are sufficiently large, the relative error of the final result between different relaxation calculations is controlled to be less than a few percent. For each finite spin-boson system, a multilayer scheme is set up where the number of the layers is chosen (empirically) 9256

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

is repeated until the contracted basis functions between different L3 single particle groups are balanced. After the multilayer hierarchy is set up, the wave function is represented in the form of (1.1) with sufficient number of configurations for each layer. In this work the converged number of configurations for each lower layer is typically 103− 104, and 104−105 for the top layer. To obtain each eigenstate, the iterative multilayer improved relaxation method discussed in section II is employed. A suitable imaginary time interval is used for each iterative relaxation step. In this work it is chosen to be longer than the inverse of the smallest energy parameter, i.e., Δ in eq 3.1. The iterations are proceeded as described in section II until convergence is reached. The size of the Krylov subspace ranges from 20 to 100 in the Lanczos/Arnoldi procedure for diagonalizing the Boltzmann operator in the top layer. At first glance this number seems surprising because only the two lowest energy eigenstates are required (in two separate calculations, i.e., not using the state-averaged approach.61) However, when quantum phase transition is studied, very high accuracy is required for the energies. A relatively large Krylov subspace ensures that any spurious contribution in the dense energy manifold is removed by the orthogonalization process of many Krylov vectors. C. Results and Discussion. We first consider a symmetric spin-boson model [ϵ = 0 in Hamiltonian (3.1)] and the nonadiabatic regime where the characteristic frequency ωc is much bigger than the two-level nonadiabatic coupling parameter Δ (we use atomic units where ℏ = 1). The scaling limit is reached when ωc/Δ → ∞, whereas for demonstration purpose we consider a large value of ωc/Δ = 40. The energy splitting between the ground and the first excited states, E1 − E0, is obtained from two separate improved relaxation calculations. As the number of modes increases, the splitting becomes smaller. This is true even with zero system-bath coupling, α = 0, due to the frequency decrease for the discrete modes at the lower end of the spectral density. For the results discussed below, a sufficient number of modes is used so that the frequency of the first mode is always smaller than the nonadiabatic coupling parameter Δ even with the smallest number of bath modes considered in this work. In the case of α = 0, the energy splitting is just that of the first mode. Thus, we adopt a scaled energy splitting, Nb(E1 − E0)/Δ, as used in the previous DMRG work,71 where Nb is the number of discrete bath modes in the calculation. For zero system-bath coupling this scaled splitting is essentially a constant versus Nb because the frequency of the first mode is roughly inversely proportional to Nb. One expects a similar dependence for a weak system-bath coupling, as shown for α = 0.05 in Figure 1a. As the coupling strength increases, the two-level subsystem participates in the energy splitting and changes this dependence. As shown in Figure 1a, the scaled energy splitting increases with the increase in Nb. For α < 0.5 the extent of this increase becomes more significant for a larger Kondo parameter α, which is due to the mixing of the two-level subsystem with the phonon bath. The trend is changed for α > 0.5, as shown in Figure 1b. In particular, when α approaches the critical value α* (slightly larger than 1) the scaled energy splitting Nb(E1 − E0)/ Δ becomes again less dependent on Nb. This signals a new scaling behavior where the scaled energy splitting is close to the concept of tunneling splitting between two local wells. Passing α*, the scaled energy splitting becomes less as Nb increases. This is shown in Figure 1c, which indicates that in the infinite bath limit the energy splitting/tunneling splitting disappears,

Figure 1. Dependence of the scaled energy splitting between the ground and the first excited states, Nb(E1 − E0)/Δ, on the number of bath modes Nb. The characteristic frequency of the bath is ωc/Δ = 40.

rendering a pair of degenerate states with each one localized on one side. In another word, in the condensed phase limit, a new quantum phase is formed where localization is present. The critical coupling α* can be estimated by examining the dependence of Nb(E1 − E0)/Δ versus Nb for different Kondo parameters. As shown in Figure 2a with Nb < 1000, for α = 1.01 the scaled splitting Nb(E1 − E0)/Δ increases with the increase in Nb, suggesting a still delocalized phase. When α increases to 1.02, Nb(E1 − E0)/Δ decreases with the increase in Nb, suggesting now a localized phase. The transition approximately occurs at α = 1.015, where the dependence of Nb(E1 − E0)/Δ on Nb is almost flat. This is the boundary for the quantum phase transition, which has been used in the previous work for the pure Ohmic (linear) spectral density (without exponential cutoff).71 The thermodynamic analysis here also agrees with the time-dependent ML-MCTDH study reported previously.36 The extracted critical coupling depends on the number of bath modes. As shown in Figure 2b, the transition boundary is closer to α = 1.01 for 1000−2000 modes. The infinite bath limit may then be extrapolated from a set of α* values with different Nb. For a multilayer improved relaxation calculation there are several convergence parameters, e.g., the number of SPFs (which determines the size of configuration space) for each layer, the number of primitive basis functions, the imaginary 9257

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

Table 2. Convergence of the Energy Splitting with Respect to the Length of Imaginary Time (τ) for a Three-Layer Improved Relaxation Calculation on a Spin-Boson Model with 500 Bath Modes imaginary time Nb(E1 − E0)/Δ

time interval for each iteration, etc. Similar to the single layer case,57 we found that the ground state is always easy to converge even with a relatively small configuration space. The excited state, however, requires that there are enough primitive basis functions and configurations to represent the Hilbert space of the wave function. Table 1 shows an example for a 500 Table 1. Convergence of the Energy Splitting with Respect to the Number of SPFs for a Three-Layer Improved Relaxation Calculation on a Spin-Boson Model with 500 Bath Modes 4 2.051

12 5 2.052

14 6 2.051

3 2.067

τΔ = 2 2.051

τΔ = 3 2.052

block form of improved relaxation using the single layer MCTDH method, 74 where both states are computed simultaneously with the same set of SPFs. This works well if the SPFs associated with the two (or more) eigenstates are similar, which is true for a weaker system-bath coupling strength (a smaller Kondo parameter α) of the spin-boson model. A good indicator for this situation is that the lowest energy eigenvalues are very close between two different multilayer improved relaxation calculations with different target states. When the Kondo parameter reaches the critical coupling strength α ∼ 1, the quantum phase transition renders that the ground and the first excited states have quite different characteristics. In this case, the calculation targeting the first excited state generates a noticeable lower eigenvalue of this (first excited) state than that obtained from the calculation targeting the ground state, whereas the ground state energy eigenvalue obtained from the former is higher than that obtained from the latter. The block form of improved relaxation will not be as efficient because it has to describe two distinctively different eigenstates via a same set of SPFs of all layers. Besides energy eigenvalues, the multilayer improved relaxation method also generates the corresponding eigenvectors in the form of (1.1). With these wave functions, many other physical quantities may be evaluated. As an example, let us evaluate the spin observables and entanglement for the ground state of an asymmetric spin-boson model [ϵ ≠ 0 in Hamiltonian (3.1)], which has been studied before using Bethe Ansatz and NRG techniques for coupling strength α < 1.72 We consider the following observables that extends the coupling strength to α > 1: the longitudinal spin magnetization ⟨σz⟩, the transverse spin magnetization ⟨σx⟩, and the von Neumann or entanglement entropy S

Figure 2. Same as Figure 1 but showing a more detailed transition once the critical coupling strength is passed.

no. SPFs top layer no. SPFs other layers Nb(E1 − E0)/Δ

τΔ = 1 2.052

4 2.052

mode spin-boson model with Kondo parameter α = 1. The imaginary time interval τ for each relaxation iteration is τΔ = 2, and the number of the primitive basis functions is converged. For this three layer improved relaxation calculation at least 12 SPFs are required for each bath single particle of the top layer, resulting in a total of 124 × 2 = 41 472 top layer configurations. This ensures that the Hamiltonian matrix in the top layer is big enough to for an accurate estimate of the eigenvalues (i.e., to resolve a couple of states within a dense energy manifold.) On the other hand, the number of the SPFs for the single particles of other layers need not be so large. Table 1 shows that 4 SPFs are sufficient to achieve convergence. If this number is not big enough, e.g., 3 SPFs, the calculation is not fully converged even with a larger number of (14) top layer SPFs. Furthermore, Table 2 shows the convergence behavior with respect to the imaginary time interval τ for each relaxation iteration, where 12 SPFs are used for each top layer single particle and 4 SPFs for each single particle of other layers. It is clear that the result is good with any reasonable choice of τ. When calculating the energy splitting, we have performed a separate calculation for each target eigenstate. There exists a

S = −p+ log 2 p+ − p− log 2 p−

(3.7a)

where 2p± = 1 ±

⟨σx⟩2 + ⟨σy⟩2 + ⟨σz⟩2

(3.7b)

and for the spin-boson model, ⟨σy⟩ = 0. Figure 3 shows ⟨σz⟩, ⟨σx⟩, and S for the ground state of the spin-boson model as a function of the Kondo parameter α, where a bath of 1000 modes with ωc/Δ = 40 is used in the calculation. The longitudinal spin magnetization ⟨σz⟩ is proportional to the electronic energy bias ϵ in the limit of ϵ → 0. Hence, the derivative d⟨σz⟩/dϵ exists and converges for small enough ϵ. Various finite ϵ’s generate different dependence of ⟨σz⟩ on α in the nonlinear response regime. In contrast, the transverse spin magnetization ⟨σx⟩ is almost identical for different ϵ’s used in the study, indicating that for this property ϵ is already small enough to give the converged value. Finally, the entanglement entropy S first increases versus α in the delocalized phase until it reaches a plateau, and then it starts to decrease with a further increase in α. 9258

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

Figure 4. Same as Figure 3 but for Nb = 10 000.

Figure 3. Longitudinal spin magnetization ⟨σz⟩, the transverse spin magnetization ⟨σx⟩, and the entanglement entropy S as a function of the Kondo parameter α for different bias ϵ. The number of bath modes is Nb = 1000. The characteristic frequency of the bath is ωc/Δ = 40.

multilayer improved relaxation will then be a valuable tool to study important problems in physics and chemistry. Still, more work needs to be done to make the approach more efficient and robust. Furthermore, in this work we adopt the improved relaxation method of Meyer and co-workers.57 Another possibility for calculating the energy eigenstates is to use the (modified) Lanczos method directly within the ML-MCTDH framework.75 To simplify the implementation, though, a state average approach is more desirable.75 This will be examined in the future in combination with treating identical particles in the second quantized framework.41

To obtain the condensed phase limit, the above analysis needs to be carried out with respect to different number of bath modes. As an illustration, Figure 4 shows the same observables as in Figure 3 but for a bath of 10000 modes. Except for ⟨σx⟩, other observables show visible differences. To obtain the infinite bath limit of these physical observables, a proper extrapolation scheme needs to be established to investigate a series of spin-boson systems with different numbers of bath modes. Sometimes appropriate order parameters (such as the scaled energy splitting used above) also need to be defined. Such a study on ⟨σz⟩, ⟨σx⟩, and S will be carried out in the future.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (575)646-3473. Notes

The authors declare no competing financial interest.

IV. CONCLUDING REMARKS Thus, the generalization of the improved relaxation approach,57 with the use of ML-MCTDH imaginary time propagation, appears to be a promising approach for obtaining energy eigenstates of large systems. When treating a structured Hamiltonian such as the spin-boson model illustrated above, the ML-MCTDH method can handle many more degrees of freedom than existing approaches such as NRG or DMRG. The



ACKNOWLEDGMENTS This work has been supported by the National Science Foundation CHE-1012479 and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. 9259

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A



Article

(25) Wang, H.; Thoss, M. Self-Consistent Hybrid Approach for Simulating Electron Transfer Reactions in Condensed Phases. Isr. J. Chem. 2002, 42, 167−182. (26) Wang, H.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory. J. Chem. Phys. 2003, 119, 1289−1299. (27) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The MultiConfigurational Time-Dependent Hartree Approach. Chem. Phys. Lett. 1990, 165, 73−78. (28) Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. Wave-Packet Dynamics within the Multiconfiguration Hartree Framework: General Aspects and Application to NOCl. J. Chem. Phys. 1992, 97, 3199− 3213. (29) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The Multiconfiguration Time-Dependent Hartree (MCTDH) Method: A Highly Efficient Algorithm for Propagating Wavepackets. Phys. Rep. 2000, 324, 1−105. (30) Meyer, H.-D.; Worth, G. A. Quantum Molecular Dynamics: Propagating Wavepackets and Density Operators Using the Multiconfiguration Time-Dependent Hartree (MCTDH) Method. Theor. Chem. Acc. 2003, 109, 251−267. (31) Thoss, M.; Kondov, I.; Wang, H. Theoretical Study of Ultrafast Heterogeneous Electron Transfer Reactions at Dye-Semiconductor Interfaces. Chem. Phys. 2004, 304, 169−181. (32) Thoss, M.; Wang, H. Quantum Dynamical Simulation of Ultrafast Molecular Processes in the Condensed Phase. Chem. Phys. 2006, 322, 210−222. (33) Wang, H.; Thoss, M. Quantum Mechanical Evaluation of the Boltzmann Operator in Correlation Functions for Large Molecular Systems: A Multilayer Multi-Configuration Time-Dependent Hartree Approach. J. Chem. Phys. 2006, 124, 034114. (34) Wang, H.; Thoss, M. Quantum Dynamical Simulation of Electron-Transfer Reactions in an Anharmonic Environment. J. Phys. Chem. A 2007, 111, 10369−10375. (35) Craig, I. R.; Thoss, M.; Wang, H. Proton Transfer Reactions in Model Condensed-Phase Environments: Accurate Quantum Dynamics Using Multilayer Multiconfiguration Time-Dependent Hartree Approach. J. Chem. Phys. 2007, 127, 144503. (36) Wang, H.; Thoss, M. From Coherent Motion to Localization: Dynamics of the Spin-Boson Model at Zero Temperature. New J. Phys. 2008, 10, 115005. (37) Wang, H.; Thoss, M. From Coherent Motion to Localization. II. Dynamics of the Spin-Boson Model with Sub-Ohmic Spectral Density at Zero Temperature. Chem. Phys. 2010, 370, 78−86. (38) Zhou, Y.; Shao, J.; Wang, H. Dynamics of Electron Transfer in Complex Glassy Environment Modeled by the Cole-Davidson Spectral Density. Mol. Phys. 2012, 110, 581−594. (39) Wang, H.; Shao, S. Dynamics of a Two-Level System Coupled to a Bath of Spins. J. Chem. Phys. 2012, 137, 22A504. (40) Frenkel, J. Wave Mechanics; Clarendon Press: Oxford, U.K., 1934. (41) Wang, H.; Thoss, M. Numerically Exact Quantum Dynamics for Indistinguishable Particles: The Multilayer Multiconfiguration TimeDependent Hartree Theory in Second Quantization Representation. J. Chem. Phys. 2009, 131, 024114. (42) Wang, H.; Pshenichnyuk, I.; Härtle, R.; Thoss, M. Numerically Exact, Time-Dependent Treatment of Vibrationally Coupled Electron Transport in Single-Molecule Junctions. J. Chem. Phys. 2011, 135, 244506. (43) Albrecht, K. F.; Wang, H.; Mühlbacher, L.; Thoss, M.; Komnik, A. Bistability Signatures in Nonequilibrium Charge Transport Through Molecular Quantum Dots. Phys. Rev. B 2012, 86, 081412. (44) Wang, H.; Thoss, M. Numerically Exact, Time-Dependent Study of Correlated Electron Transport in Model Molecular Junctions. J. Chem. Phys. 2013, 138, 134704. (45) Wang, H.; Thoss, M. Multilayer Multiconfiguration TimeDependent Hartree Study of Vibrationally Coupled Electron Transport Using the Scattering-State Representation. J. Phys. Chem. A 2013, 117, 7431.

REFERENCES

(1) Harris, R. A.; Silbey, R. On the Stabilization of Optical Isomers through Tunneling Friction. J. Chem. Phys. 1983, 78, 7330−7333. (2) Garg, A.; Onuchic, J. N.; Ambegaokar, V. Effect of Friction on Electron Transfer in Biomolecules. J. Chem. Phys. 1985, 83, 4491− 4503. (3) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, W. Dynamics of the Dissipative Two-State System. Rev. Mod. Phys. 1987, 59, 1−85. (4) 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−3375. (5) Hartmann, L.; Goychuk, I.; Grifoni, M.; Hänggi, P. Driven Tunneling Dynamics: Bloch-Redfield Theory Versus Path-Integral Approach. Phys. Rev. E 2000, 61, R4687−R4690. (6) Ankerhold, J.; Lehle, H. Low Temperature Electron Transfer in Strongly Condensed Phase. J. Chem. Phys. 2004, 120, 1436−1449. (7) Wang, H.; Song, X.; Chandler, D.; Miller, W. H. Semiclassical Study of Electronically Nonadiabatic Dynamics in the CondensedPhase: Spin-Boson Problem with Debye Spectral Density. J. Chem. Phys. 1999, 110, 4828−4840. (8) Stock, G.; Thoss, M. Classical Description of Nonadiabatic Quantum Dynamics. Adv. Chem. Phys. 2005, 131, 243−375. (9) Stier, W.; Prezhdo, O. V. Thermal Effects in the Ultrafast Photoinduced Electron Transfer from a Molecular Donor Anchored to a Semiconductor Acceptor. Isr. J. Chem. 2002, 42, 213−224. (10) Stier, W.; Duncan, W. R.; Prezhdo, O. V. Thermally Assisted Sub-10 fs Electron Transfer in Dye-Sensitized Nanocrystalline TiO2 Solar Cells. Adv. Mater. 2004, 16, 240−244. (11) Duncan, W. R.; Stier, W. M.; Prezhdo, O. V. Ab Initio Nonadiabatic Molecular Dynamics of the Ultrafast Electron Injection Across the Alizarin-TiO2 Interface. J. Am. Chem. Soc. 2005, 127, 7941−7951. (12) Sun, X.; Wang, H.; Miller, W. H. Semiclassical Theory of Electronically Nonadiabatic Dynamics: Results of a Linearized Approximation to the Initial Value Representation. J. Chem. Phys. 1998, 109, 7064−7074. (13) Martin-Fierro, E.; Pollak, E. Semiclassical Initial Value Series Solution of the Spin Boson Problem. J. Chem. Phys. 2007, 126, 164108. (14) Makri, N.; Makarov, D. E. Tensor Propagator for Iterative Quantum Time Evolution of Reduced Density Matrices. I. Theory. J. Chem. Phys. 1995, 102, 4600−4610. (15) Sim, E.; Makri, N. Path Integral Simulation of Charge Transfer Dynamics in Photosynthetic Reaction Centers. J. Phys. Chem. B 1997, 101, 5446−5458. (16) Mühlbacher, L.; Egger, R. Crossover from Nonadiabatic to Adiabatic Electron Transfer Reactions: Multilevel Blocking Monte Carlo Simulations. J. Chem. Phys. 2003, 118, 179−191. (17) Egger, R.; Weiss, U. Quantum Monte Carlo Simulation of the Dynamics of the Spin-Boson Model. Z. Phys. B 1992, 89, 97−107. (18) Egger, R.; Mak, C. H. Low-Temperature Dynamical Simulation of Spin-Boson Systems. Phys. Rev. B 1994, 50, 15210−15220. (19) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (20) Bulla, R.; Tong, N.-G.; Vojta, M. Numerical Renormalization Group for Bosonic Systems and Application to the Sub-Ohmic SpinBoson Model. Phys. Rev. Lett. 2003, 91, 170601. (21) Anders, F. B.; Schiller, A. Spin Precession and Real-Time Dynamics in the Kondo Model: Time-Dependent Numerical Renormalization-Group Study. Phys. Rev. B 2006, 74, 245113. (22) Anders, F. B.; Bulla, R.; Vojta, M. Equilibrium and Nonequilibrium Dynamics of the Sub-Ohmic Spin-Boson Model. Phys. Rev. Lett. 2007, 98, 210402. (23) Wang, H.; Thoss, M.; Miller, W. H. Systematic Convergence in the Dynamical Hybrid Approach for Complex Systems: A Numerically Exact Methodology. J. Chem. Phys. 2001, 115, 2979−2990. (24) Thoss, M.; Wang, H.; Miller, W. H. Self-Consistent Hybrid Approach for Complex Systems: Application to the Spin-boson Model with Debye Spectral Density. J. Chem. Phys. 2001, 115, 2991−3005. 9260

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261

The Journal of Physical Chemistry A

Article

(46) Wilner, E. Y.; Wang, H.; Cohen, G.; Thoss, M.; Rabani, E. Bistability in a Nonequilibrium Quantum System with ElectronPhonon Interactions. Phys. Rev. B 2013, 88, 045137. (47) Kondov, I.; Cizek, M.; Benesch, C.; Wang, H.; Thoss, M. Quantum Dynamics of Photoinduced Electron-Transfer Reactions in Dye-Semiconductor Systems: First-Principles Description and Application to Coumarin 343-TiO2. J. Phys. Chem. C 2007, 111, 11970− 11981. (48) Thoss, M.; Kondov, I.; Wang, H. Correlated Electron-Nuclear Dynamics in Ultrafast Photoinduced Electron-Transfer Reactions at Dye-Semiconductor Interfaces. Phys. Rev. B 2007, 76, 153313. (49) Vendrell, O.; Meyer, H.-D. Multilayer Multiconfiguration Timedependent Hartree method: Implementation and Applications to a Henon-Heiles Hamiltonian and to Pyrazine the Description of Polyatomic Molecules. J. Chem. Phys. 2011, 134, 044135. (50) Kondov, I.; Wang, H.; Thoss, M. Theoretical Study of Ultrafast Heterogeneous Electron Transfer Reactions at Dye-Semiconductor Interfaces: Coumarin 343 at Titanium Oxide. J. Phys. Chem. A 2006, 110, 1364−1374. (51) Wang, H.; Thoss, M. Nonperturbative Quantum Simulation of Time-Resolved Nonlinear Spectra: Methodology and Application to Electron Transfer Reactions in the Condensed Phase. Chem. Phys. 2008, 347, 139−151. (52) Egorova, D.; Gelin, M. F.; Thoss, M.; Wang, H.; Domcke, W. Effects of Intense Femtosecond Pumping on Ultrafast ElectronicVibrational Dynamics in Molecular Systems with Relaxation. J. Chem. Phys. 2008, 129, 214303. (53) Velizhanin, K. A.; Wang, H.; Thoss, M. Heat Transport Through Model Molecular Junctions: A Multilayer Multiconfiguration Time-Dependent Hartree Approach. Chem. Phys. Lett. 2008, 460, 325−330. (54) Velizhanin, K. A.; Wang, H. Dynamics of Electron Transfer Reactions in the Presence of Mode-Mixing: Comparison of a Generalized Master Qquation Approach with the Numerically Exact Simulation. J. Chem. Phys. 2009, 131, 094109. (55) Manthe, U. A Multilayer Multiconfigurational Time-Dependent Hartree Approach for Quantum Dynamics on General Potential Energy Surfaces. J. Chem. Phys. 2008, 128, 164116. (56) Manthe, U. Layered Discrete Variable Representations and Their Application Within the Multiconfigurational Time-Dependent Hartree Approach. J. Chem. Phys. 2009, 130, 054109. (57) Meyer, H.-D.; Quéré, F. L.; Léonard, C.; Gatty, F. Calculation and Selective Population of Vibrational Levels with the Multiconfiguration Time-Dependent Hartree (MCTDH) Algorithm. Chem. Phys. 2006, 329, 179−192. (58) Manthe, U.; Matzkies, F. Iterative Diagonalization within the Multi-Configurational Time-Dependent Hartree Approach: Calculation of Vibrationally Excited States and Reaction Rates. Chem. Phys. Lett. 1996, 252, 71−76. (59) Matzkies, F.; Manthe, U. A Multi-Configurational TimeDependent Hartree Approach to the Direct Calculation of Thermal Rate Constants. J. Chem. Phys. 1997, 106, 2646−2653. (60) Hinze, J. MC-SCF. I. The Multiconfiguration Self-Consistent Field Method. J. Chem. Phys. 1973, 73, 6424−6432. (61) Manthe, U. The State Averaged Multiconfigurational Timedependent Hartree Approach: Vibrational State and Reaction Rate Calculations. J. Chem. Phys. 2008, 128, 064108. (62) Young, D. M.; Mai, T.-Z. In Iterative Methods for Large Linear Systems; Kincaid, D. R., Haynes, L. J., Eds.; Academic Press: San Diego, 1990. (63) Weiss, U. Quantum Dissipative Systems, 2nd ed.; World Scientific: Singapore, 1999. (64) Marcus, R. A.; Sutin, N. Electron Transfers in Chemistry and Biology. Biochim. Biophys. Acta 1985, 811, 265−322. (65) Suárez, A.; Silbey, R. Properties of a Macroscopic System as a Thermal Bath. J. Chem. Phys. 1991, 95, 9115−9121. (66) Weiss, U.; Grabert, H.; Linkwitz, S. Influence of Friction and Temperature on Coherent Quantum Tunneling. J. Low Temp. Phys. 1987, 68, 213−244.

(67) Spohn, H. Ground State(s) of the Spin-Boson Hamiltonian. Commun. Math. Phys. 1989, 123, 277−304. (68) Bray, A.; Moore, N. Influence of Dissipation on Quantum Coherence. Phys. Rev. Lett. 1982, 49, 1545−1549. (69) Chakravarty, S. Quantum Fluctuation in the Tunneling between Superconductors. Phys. Rev. Lett. 1982, 49, 681−684. (70) Kehrein, S. K.; Mielke, A.; Neu, P. Flow Qquations for the SpinBoson Problem. Z. Phys. B 1996, 99, 269−280. (71) Wong, H.; Chen, Z.-D. Density Matrix Renormalization Group Approach to the Spin-Boson Model. Phys. Rev. B 2008, 77, 174305. (72) Hur, K. L. Entanglement Entropy, Decoherence, and Quantum Phase Transitions of a Dissipative Two-Level System. Ann. Phys. 2008, 323, 2208−2240. (73) Wang, H.; Thoss, M. Theoretical Study of Ultrafast Photoinduced Electron Transfer Processes in Mixed-Valence Systems. J. Phys. Chem. A 2003, 107, 2126−2136. (74) Doriol, L. J.; Gatti, F.; Iung, C.; Meyer, H.-D. Computation of Vibrational Energy Levels and Eigenstates of Fluroform Using the Multiconfiguration Time-Dependent Hartree Method. J. Chem. Phys. 2008, 129, 224109. (75) Hammer, T.; Manthe, U. Intramolecular Proton Transfer in Malonaldehyde: Accurate Multilayer Multiconfiguration Time-Dependent Hartree Calculations. J. Chem. Phys. 2011, 134, 224305.

9261

dx.doi.org/10.1021/jp503351t | J. Phys. Chem. A 2014, 118, 9253−9261