Quantum Phase Transition in the Spin-Boson Model: A Multilayer

Subscriber access provided by TULANE UNIVERSITY

Article

Quantum Phase Transition in the Spin-Boson Model: A Multilayer Multiconfiguration Time-Dependent Hartree Study Haobin Wang, and Jiushu Shao J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b11136 • Publication Date (Web): 08 Feb 2019 Downloaded from http://pubs.acs.org on February 10, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Quantum Phase Transition in the Spin-Boson Model: A Multilayer Multiconfiguration Time-Dependent Hartree Study Haobin Wang∗ Department of Chemistry, University of Colorado Denver, Denver, CO 80217-3364, USA, and Beijing Computational Science Research Center, No.10 East Xibeiwang Road, Haidian District, Beijing 100193, China Jiushu Shao College of Chemistry and Center for Advanced Quantum Studies, Beijing Normal University; Key Laboratory of Theoretical and Computational Photochemistry, Ministry of Education, Beijing Normal University, Beijing 100875, China

∗

Email: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 41

ABSTRACT

The multilayer improved relaxation is applied to study the delocalization-localization transition in the spin-boson model at zero temperature – a well-known example of quantum phase transition. Calculations of energy eigenstates are obtained by iteratively diagonalizing the matrix of the Boltzmann operator in the top layer representation, using a Lanczos/Arnoldi method while relaxing the single particle functions of all layers using the multilayer multiconfiguration time-dependent Hartree imaginary time propagation. Two properties are used to examine the quantum phase transition: the energy splitting for the lowest pair of eigenstates and the magnetic susceptibility. Consistent findings are obtained with appropriate scaling parameters.

2


Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I.

INTRODUCTION

The description of many-body quantum effects for processes in the condensed phase is essential to the understanding of many important phenomena. Physical models with a few parameters are particularly useful in such studies due to their universality and compactness. The spin-boson Hamiltonian is such an example that can be used to model diversified phenomena in electron transfer reactions,1 hydrogen tunneling,2 macroscopic quantum coherence,3 quantum dissipation,4,5 , quantum optics, and quantum computation. When treating such physical models, quantum perturbation theories may be used to describe the reduced dynamics of a system with weak interaction.4,6–10 Another important direction is the development of numerically exact approaches, e.g. numerical path integral approaches using Feynman-Vernon influence functional,11–15 the numerical renormalization group (NRG) theory,16–18 the stochastic field method,19–25 and the multilayer multiconfiguration timedependent Hartree (ML-MCTDH) theory.26–29 In ML-MCTDH the wave function is expressed by an iterative (or recursive), layered expansion, |Ψ(t)i =

XX j1

...

j2

X

Aj1 j2 ...jp (t)

p Y

κ=1

jp

(κ)

|ϕjκ (t)i,

(1.1a)

Q(κ) (κ) |ϕjκ (t)i

=

XX i2

i1

...

X

κ Biκ,j (t) 1 i2 ...iQ(κ)

Y q=1

iQ(κ)

(κ,q)

|viq

(t)i,

(1.1b)

M (κ,q) (κ,q) |viq (t)i

=

XX α1

...

α2

X

αM (κ,q)

q Cακ,q,i (t) 1 α2 ...αM (κ,q)

Y

γ=1

...

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

(1.1c)

κ,q,i

κ where Aj1 j2 ...jp (t), Biκ,j (t), Cα1 α2q...αM (κ,q) (t) and so on are the expansion coefficients for 1 i2 ...iQ(κ)

(κ)

(κ,q)

the first, second, third, ..., layers, respectively; |ϕjκ (t)i, |viq

(κ,q,γ)

(t)i, |ξαγ

(t)i, ..., are the

“single 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 time-independent configurations, each of which may be defined by several Cartesian degrees of freedom. The variational parameters within the ML-MCTDH theoretical framework are dynamically optimized through the use of Dirac-Frenkel variational principle,30 ∂ ˆ hδΨ(t)|i ∂t − H|Ψ(t)i = 0, which results in a set of coupled, nonlinear differential equations

for the expansion coefficients of all layers.26–29 The flexibility of the ML-MCTDH wave function provides a great advantage in the variational functional, which results in a tremendous 3



Page 4 of 41

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.31–54 The ML-MCTDH has also been implemented recursively by Manthe.55,56 The applicability of the ML-MCTDH theory is not limited to simulating real time quantum dynamics, as evidenced by applications of MCTDH for calculating eigenstates of the Hamiltonian57–59 or the thermal flux operator.59,60 Along this line the multilayer improved relaxation method was developed.51 Although for large systems the calculation of many eigenstates is less interesting due to a large density of states, it may be useful to calculate a few lowest energy eigenstates under critical conditions so as to predict important thermodynamic properties of the system. Combined with appropriate scaling parameters, such a study may be useful for extracting limits in the condensed phase. In this paper we discuss the application of the multilayer improved relaxation to the study of quantum phase transition in the spin-boson model. Section II summarizes the theory and practical implementation of the method as well as the physical model under investigation. Section III discusses our findings. Section IV concludes.

II.

THEORY AND MODEL A.

Multilayer improved relaxation

To compute energy eigenstates with the multilayer wave function ansatz, Eq. (1.1), one may perform the time-independent variation

ˆ δ hΨ|H|Ψi −E(hΨ|Ψi − 1) −

p X X κ=1 j,l

−

p Q(κ) X XX

−

p Q(κ) M (κ,q) X X X X

κ=1 q=1

κ=1 q=1

−...

(κ,q)

(κ,q)

ǫ˜j,l (hvj

j,l

γ=1

j,l

(κ)

(κ,q)

|vl

(κ)

i − δij )

˜˜ǫ(κ,q,γ) (hξj(κ,q,γ)|ξ (κ,q,γ)i − δij ) j,l l

= 0,

4

(κ)

ǫj,l (hϕj |ϕl i − δij )


(2.1)

Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(κ)

(κ,q)

where E is the Lagrange multiplier to keep the overall wave function normalized, ǫj,l , ˜ǫj,l , (κ,q,γ) ǫ˜˜j,l , and so on are Lagrange multipliers to keep the SPFs of the first, second, third, and

lower layers orthonormal. Variation in the top layer gives the standard eigen-equation HJL AL ≡

X L

ˆ L iAL = EAJ , hΦJ |H|Φ

(2.2)

where we have denoted a configuration |ΦJ i as the Hartree product |ΦJ i ≡

(κ) κ=1 |ϕjκ i

Qp

and used the collective index, Aj1 j2 ...jp ≡ AJ . Solving this eigen-equation thus gives the P eigen-energy E and the eigenvector |Ψi = J AJ |ΦJ i in terms of the top layer expansion P coefficients AJ and the SPFs of all the layers (since |Ψi = J AJ |ΦJ i depends on them).

Variation (2.1) with respect to the SPFs of all the layers, after eliminating all the Lagrange

multipliers, results in the same form ˆ (κ) |ϕ(κ) i = 0, [1 − Pˆ (κ) ]hHi

(2.3a)

(κ,q) ˆ (κ,q) (κ,q) i = 0, [1 − PˆL2 ]hHi L2 |v

(2.3b)

(κ,q,γ) ˆ (κ,q,γ) (κ,q,γ) [1 − PˆL3 ]hHi |ξ i = 0, L3

(2.3c)

... (κ,q)

(κ,q,γ)

ˆ (κ) , hHi ˆ ˆ where hHi L2 , hHiL3

, and so on are mean-field operators for the first, second, (κ,q) (κ,q,γ) third, and lower layers, respectively; Pˆ (κ) , PˆL2 , PˆL3 , and so on are projection operators in the single particle space using the SPFs of the corresponding layers; and |ϕ(κ) i = (κ)

(κ)

(κ,q)

{|ϕ1 i, |ϕ2 i, ...}T , |v(κ,q) i = {|v1

(κ,q)

i, |v2

(κ,q,γ)

i, ...}T , |ξ (κ,q,γ) i = {|ξ1

(κ,q,γ)

i, |ξ2

i, ...}T , and

so on 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 has been detailed in the previous ML-MCTDH publications.26–29 The above set of equations, Eqs.(2.2)-(2.3), are highly nonlinear and very difficult to converge with iteration. To get around this difficulty the MCTDH “improved relaxation” method58 was generalized to the ML-MCTDH framework.51 Since the stationary solution in (2.3) remains valid under any linear transformation, one may insert the pseudo-inverse (κ,q)

(κ,q,γ)

of the reduced density matrices ρˆ(κ) , ̺ˆL2 , ̺ˆL3

, ..., to the equation of the corresponding

layer ˆ (κ) |ϕ(κ) i = 0, [1 − Pˆ (κ) ][ˆ ρ(κ) ]−1 hHi 5


(2.4a)


Page 6 of 41

(κ,q) (κ,q) ˆ (κ,q) |v(κ,q) i = 0, [1 − PˆL2 ][ˆ ̺L2 ]−1 hHi L2

(2.4b)

(κ,q,γ) (κ,q,γ) ˆ (κ,q,γ) |ξ (κ,q,γ) i = 0, [1 − PˆL3 ][ˆ ̺L3 ]−1 hHi L3

(2.4c)

... and realize that these can be related to the stationary solution of the following evolution of the SPFs with imaginary time, i.e. ˆ (κ) |ϕ(κ) i, |ϕ˙ k (τ )iL2 coefficients = −[1 − Pˆ (κ) ][ˆ ρ(κ) ]−1 hHi

(2.5a)

(κ,q) (κ,q) ˆ (κ,q) |v(κ,q) i, ̺L2 ]−1 hHi |v˙ (κ,q) (τ )iL3 coefficients = −[1 − PˆL2 ][ˆ L2

(2.5b)

(κ,q,γ)

|ξ˙

(κ,q,γ) (κ,q,γ) ˆ (κ,q,γ)|ξ (κ,q,γ) i, (τ )iL4 coefficients = −[1 − PˆL3 ][ˆ ̺L3 ]−1 hHi L3

(2.5c)

... In the notation used in Eq. (2.5) 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 taken for the κ coefficient function Biκ,j (t) defined in (1.1b), and so on. Finally, the imaginary time is 1 i2 ...iQ(κ)

defined as τ = −it. Thus, the approximate solution to Eq. (2.3) is to use the ML-MCTDH method to propagate the SPFs in imaginary time (relaxation) by applying Eq. (2.5) until the time-derivatives are sufficiently small. This is then coupled to the solution of the top layer eigen-equation (2.2) in an iterative way. Overall, the multilayer improved relaxation algorithm is realized as follows.51 An initial guess is first generated, e.g., by diagonalizing a zeroth-order Hamilˆ 0 and choosing an appropriate target (ground or excited) state. Then one proceeds tonian H to the iterations of multilayer improved relaxation. In each iteration step, the Hamiltonian matrix of the top layer is first constructed based on the current SPFs. Then this Hermitian matrix is diagonalized as shown in Eq. (2.2). In our work an Arnoldi/Lanczos method on the Boltzmann operator is employed, in spirit similar to the work of Manthe59,61 but only applying the Krylov iteration to the linear part of the top layer A-coefficients. One may also use Davidson’s method as in the previous work by Meyer and coworkers58 with a preconditioner, e.g., Jacobi (diagonal), Gauss-Seidel, or symmetric successive over relaxation (SSOR)62 to accelerate the Krylov iteration. After diagonalizing the top layer Hamiltonian matrix, the A-coefficients of the target eigenstate is fixed while the SPFs of all layers are 6


Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 is achieved for the target state of interest. In this work, numerical convergence criterion is set by requiring that both the energy and the reduced density agree within 10−8 in two consecutive iterations.

B.

The Spin-Boson Model and Wilson’s Logarithmic Discretization

The spin-boson Hamiltonian,4,5 in the context of electron transfer theory, comprises two electronic states, |φ1 i and |φ2 i, linearly coupled to a bath of harmonic oscillators. In massweighted coordinates the Hamiltonian reads X ∆ 1X 2 ǫ (pi + ωi2 qi2 ) + σz ci qi , H = σz + σx + 2 2 2 i i

(2.6)

where σx and σz are Pauli matrices σx = |φ1 ihφ2 | + |φ2ihφ1 |,

(2.7a)

σz = |φ1 ihφ1 | − |φ2 ihφ2 |.

(2.7b)

The properties of the bath that influence the dynamics of the two-state subsystem are specified by the spectral density function4,5 J(ω) =

π X c2j δ(ω − ωj ). 2 j ωj

(2.8)

In this work we use the Ohmic spectral density that is often employed in NRG studies   π αω (0 < ω ≤ ωc ) 2 J(ω) = (2.9) 0 (ω > ωc ), where α is the dimensionless Kondo parameter that characterizes the system-bath coupling

strength and ωc is the cutoff frequency of the bath. The relation of the spin-boson model to the Kondo problem has been discussed previously.4,5 The continuous bath spectral density of Eq. (2.9) can be discretized to the form of Eq. (2.8) via the relation26,63–65 c2j =

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


(2.10)


Normally the density of frequencies ρ(ω) is defined from the integral relation Z ωj dω ρ(ω) = j, j = 1, ..., N.

Page 8 of 41

(2.11)

0

In this work, we employ Wilson’s logarithmic discretization that is used in the NRG work. The frequencies are given by ωN

≡ ωc , ωc = , Λ

ωN −1 ... ωN −j

=

ωc , Λj

=

ωc , ΛN −1

... ω1

(2.12)

where Λ → 1+ and N is the number of discrete bath modes. It is easy to show that this discretization is related to the following density of frequencies ρ(ω) =

1 , ωlnΛ

with a slightly modified version of (2.11) Z ωc dω ρ(ω) = N − j,

(2.13)

j = 1, ..., N.

(2.14)

ωj

In some applications one chooses a variant to frequencies ωk in Eq. (2.12), i.e. using the average frequencies ω ¯k

where

Rk

dω ρ(ω)ω , ω ¯k = R k dω ρ(ω) Z

k

dω ≡

Z

(2.15a)

ωk

dω.

(2.15b)

ωk−1

Moreover, the density of frequencies ρ(ω) does not need to be a continuous function in the whole interval as Eq. (2.13). Instead, in some applications ρ(ω) is chosen to be piecewise continuous in each of the [ωk−1 , ωk ] interval, where ωk is given in (2.12). Within each such interval one is free to choose a particular ρ(ω) such that Z k dω ρ(ω) = 1.

(2.16)

For example, one may require that ρ(ω) ∝ J(ω) in each interval as done in some NRG applications. Our investigation shows that these different variants make negligible difference in our applications when the number of discrete modes is large enough. 8


Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


C.

Properties to Quantify Phase Transition

In previous studies the energy splitting between the ground and the first excited state, E1 − E0 , was often used to categorize the delocalized-localized phase transition in the spinboson model. When E1 − E0 is finite, the system possesses the parity symmetry and its ground state is delocalized. By contrast, when the ground and the first excited states become degenerate, i.e, E1 − E0 = 0, the parity symmetry is broken and its ground states are localized. This delocalized-localized or quantum phase transition occurs beyond a critical coupling strength αc .66,67 In the nonadiabatic limit (∆/ωc ≪ 1) and with a zero energy bias/level asymmetry (ǫ = 0), the phonon-dressed tunneling splitting, ∆eff , is an important parameter as compared with the energy difference between the two adjacent states for the lowest-frequency mode, i.e. ω1 in (2.12). (Here we use atomic units where h ¯ = 1.) In the weak coupling limit, ∆eff > ω1 when the number of modes N is greater than a threshold value Nm . As a result, the energy splitting E1 − E0 for the whole spin-boson system, with N finite modes, approaches the energy gap of the lowest-frequency mode E1 − E0 → ω1 =

ωc , ΛN −1

(2.17)

as N → ∞. This suggests the scaled energy splitting Es = ΛN (E1 − E0 ) as an appropriate physical variable and the existence of a “flow diagram”, as shown in Figure 1a, that is often applied in NRG studies.16 The results in Figure 1 were obtained using the multilayer improved relaxation method for ωc /∆ = 10 and the scaling factor Λ = 1.05. For zero coupling (α = 0) the scaled energy splitting is given by Es /∆ ≡ ΛN (E1 − E0 )/∆ = Λωc /∆ = 10Λ

(2.18)

for N > Nm (in this case Nm = 1 + ln10/ln1.05 = 48), as evidenced from expression (2.17) and the condition ω1 < ∆. As shown in Figure 1a, this limit is also reached in the weak coupling regime, with Nm increasing versus α. In the previous NRG study, a “crossover” scale T ∗ was defined as16 T ∗ = const. × Λ−N , ∗

(2.19)

where N ∗ is defined as the value of N where the scaled energy splitting Es = ΛN (E1 − E0 ) reaches a particular value that is less than the plateau Λωc . Ref. 16 asserts that this value 9



Page 10 of 41

is arbitrary based on the assumption that its change can be absorbed in a change of the prefactor in the above expression. Then, by virtue of renormalization group arguments the critical coupling αc can be extracted from the following expression for the dependence of T ∗ on α and ∆/ωc ∗

T ∝

∆ ωc

1/(αc −α)

,

(2.20a)

or 1 lnT = C + ln αc − α ∗

∆ ωc

.

(2.20b)

Another possibility is to examine the qualitative change in Es = ΛN (E1 − E0 ) versus N for different values of the Kondo parameter α. This has been done in a previous study employing the density matrix renormalization group (DMRG) approach,68 as well as the illustration of the multilayer relaxation method.51 The difference is the definition of the scaled energy is N(E1 − E0 ) due to the equi-distance discretization of the bath, instead

of ΛN (E1 − E0 ) for the logarithmic discretization used in the NRG study and this work. Under this approach, when α approaches the critical value αc the scaled energy splitting Es

becomes less dependent on N. This signals a new scaling behavior where the energy splitting is close to zero as the tunneling between the two localized states stops. Before and after αc , the dependence of the scaled energy splitting Es on the Kondo parameter α, distinguished by the critical value αc , displays a qualitative difference as shown in Figure 1b: for α = 1 the scaled splitting Es = ΛN (E1 − E0 ) increases with the increase in N, suggesting a delocalized phase. When α increases beyond 1.03, Es decreases with the increase in N, which indicates that in the infinite bath limit the energy splitting/tunneling splitting disappears, 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 transition approximately occurs at α = 1.02 − 1.03, where the dependence of Es on N is relatively flat. This is the boundary for the quantum phase transition, the behavior of which (as well as the quantitative values) will be discussed in more detail in the next section. It should be noted that passing the phase boundary the scaled splitting Es = ΛN (E1 − E0 ) is smaller than the value before the boundary. In particular, since it is a decreasing function of N, its value should be less than one which is the value at N = 0. To characterize the quantum phase transition one may consider a second quantity, namely, the magnetic susceptibility. By definition, it is the derivative of the longitudinal spin mag10


Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


netization hσz i for the many-body ground state |Ψ0 i with respect to the external field represented by the detuning/energy bias, dhσz i d . − = − hΨ |σ |Ψ i 0 z 0 dǫ ǫ=0 dǫ ǫ=0

(2.21)

dhσz i Again, a proper scaling versus N and ∆ is introduced, hereby as −Λ−N d(ǫ/∆) . In the ǫ=0

absence of the bath (N = 0) the scaled derivative is 1. With the exponential scaling

factor Λ−N the scaled derivative decreases exponentially versus N at zero coupling. As illustrated in Figure 2a, the remnant of this exponential decrease can be found in the weak coupling regime. As coupling becomes stronger, the decrease in the scaled derivative becomes slower. Meanwhile some curvature appears on a logarithmic scale. Quantum phase transition dhσz i changes from decreasing to increasing versus N, as evidenced occurs when −Λ−N d(ǫ/∆) ǫ=0 in Figure 2b. Just as the tunneling splitting disappears for the localized phase, − dhσz i dǫ

ǫ=0

becomes infinitely large for the localized phase.

Thus, for the example considered in this paper we examine the dependence of a scaled physical quantity on the number of the bath modes N. This way a thermodynamic limit may be extracted in the limit of an infinite bath via some proper scaling procedure. The multilayer relaxation method allows us to use a relatively small scaling parameter Λ → 1+ that is formally required by the Wilson logarithmic discretization. The results presented below employ Λ = 1.05 and 1.01. The number of bath modes in the simulation is up to 200 for the former and up to 1000 for the latter.

III.

RESULTS AND DISCUSSION

We first consider using Eq. (2.20) to extrapolate the critical coupling αc . Under perfect conditions, it is assumed that one may use an arbitrary Es /∆ to define the crossover scale in (2.19) and then (2.20). In reality, though, this cutoff value Es /∆ should not be too small because the resulting N ∗ in (2.19) would be far away from the condensed phase limit. It is also obvious that the upper limit of this cutoff value is ωc /∆. Figure 3 illustrates the extrapolation of αc that corresponds to the results displayed in Figure 1. In Figure 3a the fitting range is α ∈ [0.3, 0.8] and the cutoff value is Es /∆ = 6. A reasonably good fit to (2.20) gives αc = 1.063. This value seems correct and is also consistent with the NRG extrapolated value.16 There are, however, some uncertainties associated with 11



Page 12 of 41

different choices of Es /∆, a parameter for fitting (2.20) via definition (2.19). If Es /∆ is selected in the range of [2.5, 7.5], the fitted αc varies from 1.207 to 1.049. The fit deteriorates for either too small or too large Es . Figure 3b displays the situation for Es /∆ = 7.5 where the dependence of T ∗ on α is clearly less accurately described by (2.20). Further complications arises if the range of α, another parameter for fitting (2.20), moves to the stronger coupling regime. Figure 4 illustrates the dependence of T ∗ on α for the fitting range α ∈ [0.5, 0.85]. When Es /∆ = 5 is used, the fitted critical coupling is αc = 1.080. When Es /∆ = 2.5 is used, however, the fitted critical coupling is αc = 1.159. As discussed below, this range of uncertainty is even bigger than the change of αc versus ωc . As described in the previous section, another way to determine αc is to examine the qualitative change in Es = ΛN (E1 − E0 ) versus N for different Kondo parameters α near αc ≃ 1, as has been done earlier.68 It should be pointed out again that in DMRG work

the definition of the scaled energy is N(E1 − E0 ) because of the equi-distance discretization of the bath. When α approaches the critical value αc the scaled energy splitting becomes closer to the tunneling splitting between the two localized potential wells. Before reaching αc , the whole system has a delocalized ground state and a phonon-dressed tunneling splitting that decreases slower than the decrease in the phonon’s lowest frequency, Λ−N . Beyond αc , the phonon-dressed tunneling splitting would disappear, which is reflected in the calculation as decreasing faster than Λ−N . This signals a pair of degenerate localized states or a delocalization-localization (quantum-classical) phase transition. As shown in Figure 1b, the transition approximately occurs at α = 1.02−1.03, where the dependence of ΛN (E1 −E0 ) on N is relatively flat. A more detailed illustration is given in Figure 5 where several Es versus N plots are included. In the scaling limit ωc → ∞, a sharp transition is expected at αc = 1+ . However, for a finite ωc this transition is softer and spans a range. Based on the changes illustrated in Figure 5 and an approximate extrapolation to N → ∞, the phase transition is estimated to be around αc = 1.03. This value is smaller than that from fitting (2.20). Similar finding was reported in the previous DMRG work with an equidistance discretized bath.68 In principle the scaling factor Λ in the logarithmic discretization should be close to unity to ensure convergence. The results described so far were obtained with Λ = 1.05. This is good enough for the purpose of estimating αc with the current cutoff frequency ωc /∆ = 10. As a demonstration, Figures 6 and 7 show the results obtained with Λ = 1.01. In Fig. 6 12


Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fittings to (2.20) are illustrated for two ranges of α. Considering the uncertainties in such fitting, the αc ’s obtained for Λ = 1.01 agree with those in Figs. 3 and 4 for Λ = 1.05. With respect to the qualitative change in Es = ΛN (E1 − E0 ) versus N, Fig. 7 shows a very similar transition to that depicted in Fig. 5. The predicted αc is also around 1.03, same as obtained from Fig. 5. We also performed the calculations with a somewhat larger scaling factor Λ = 2, which was used in previous NRG studies. Figure 8 illustrates the results using two ways to extrapolate the critical coupling strength αc . In Fig. 8a fittings to (2.20) is used. The fitted αc = 1.067 is consistent with the values obtained with smaller scaling factors, as shown above, and is certainly within the error range of such a fitting protocol. Similarly, the qualitative change in Es = ΛN (E1 − E0 ) versus N also predicts αc = 1.0275 − 1.03, again the same as that obtained above with smaller scaling factors. On one hand, this indicates that a scaling of Λ = 2 is already good enough for the purpose of estimating αc , certainly within the error range of eigenstate calculations. On the other hand, it raises questions of extrapolating αc as a function of the scaling factor Λ that has been used in previous NRG studies, because our results indicate αc is relatively independent of Λ when the latter is small enough. Some other observables were proposed for studying phase transitions in the spin-boson model using Bethe ansatz and NRG techniques.69 Examples include the longitudinal spin magnetization hσz i, the transverse spin magnetization hσx i, and the von Neumann or entanglement entropy S S = −p+ log2 p+ − p− log2 p− ,

(3.1a)

where 2p± = 1 ±

q

hσx i2 + hσy i2 + hσz i2 .

(3.1b)

However, we did not find a reliable numerical criterion to use any of the above for classifying the phase transition since these quantities do not exhibit qualitative changes before and after the phase transition. Instead, since the ground-state longitudinal spin magnetization hσz i is proportional to the electronic energy bias ǫ in the limit of weak coupling and ǫ → 0, the derivative dhσz i/dǫ exists and converges for small enough ǫ. A properly defined quantity can thus be used to quantify the phase transition. This is the scaled derivative −Λ−N dhσz i d(ǫ/∆) ǫ=0

at the ground state, as shown in Fig. 2. The derivative is obtained from numerical finite

13



Page 14 of 41

difference hσz iǫ − hσz i0 hσz iǫ dhσz i ≃ = (3.2) . ǫ→0 d(ǫ/∆) ǫ=0 ǫ/∆ ǫ/∆ ǫ→0 dhσz i Comparing Fig. 1b and Fig. 2b, it is found that −Λ−N d(ǫ/∆) provides similar information ǫ=0

as the scaled energy splitting Es , only with an opposite trend: as α goes beyond αc , a tiny

bias ǫ will induce a big change in hσz i that scales more quickly than ΛN . dhσz i because it only requires calculating the ground It may seem easier to use −Λ−N d(ǫ/∆) ǫ=0 state as compared with calculating the ground and the first excited states when evaluating the energy splitting. However, it is more challenging to evaluate the numerical derivative via (3.2). Although hσz i is proportional to the electronic energy bias ǫ in the limit of ǫ → 0, various finite ǫ’s generate different dependence of hσz i on α in the nonlinear response regime. As illustrated in Fig. 9, the numerical derivative is not reliable at larger N for an ǫ too large, e.g., ǫ/∆ = 10−4. On the other hand, small numerical noise enters for a too small ǫ, e.g., ǫ/∆ = 10−7 . In most situations one needs to contemplate between these two limits and this creates numerical uncertainties. This makes the scaled energy splitting Es = ΛN (E1 − E0 ) a more robust numerical measure for the phase transition, although the two provide the same information as shown by comparing Figures 1 and 2. As discussed above, however, the method we employ to extrapolate αc , using the qualitative change in the scaled energy splitting Es = ΛN (E1 − E0 ), yields different answers from Eq. (2.20) that is used in the previous NRG extrapolation approach, even considering the numerical uncertainties. The difference of estimating the critical coupling strength αc is summarized in Table 1. Figures 10 and 11 illustrate such differences for a lower cutoff frequency ωc /∆ = 5 and a higher one ωc /∆ = 20. The fit to (2.20) generates very similar results for these two frequencies: αc = 1.085 for the former and αc = 1.076 for the latter. Such a closeness of the extrapolated αc ’s for different cutoff frequencies ωc , using the NRG scheme (2.20), were also found in the previous NRG studies.16 This is in contrast to our current approach of extrapolating αc by examining the qualitative difference of Es = ΛN (E1 − E0 ) around the critical coupling. As illustrated in Figs. 10b and 11b, the quantum phase transition occurs around αc = 1.08 − 1.10 for ωc /∆ = 5 and αc = 1.005 − 1.01 for ωc /∆ = 20. Combined with αc = 1.02 − 1.03 for ωc /∆ = 10

shown in Figs. 7, one may quantitatively predict (extrapolate) αc → 1+ in the scaling limit

ωc → ∞. This is even more evident in Figure 12 for ωc /∆ = 100, where the extrapolated 14


Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


αc = 1.0004 − 1.0005 in Figure 12b becomes very close to unity which is the scaling limit. TABLE 1: Critical coupling αc for the various parameters and the two ways of extracting it. Scaled energy splitting criterion ωc /∆

Fit to Eq. (2.20)

lower bound

upper bound

lower bound upper bound

5

1.08

1.10

1.08

1.13

10

1.02

1.03

1.08

1.12

20

1.005

1.01

1.08

1.12

100

1.0004

1.0005

1.08

1.13

On the contrary, Figure 12a illustrates the uncertainty of using Eq. (2.20) to extrapolate αc with ωc /∆ = 100. For α ∈ [0.35, 0.75] and Es /∆ = 13, a fit to (2.20) yields a value of αc = 1.13. On the other hand, for α ∈ [0.25, 0.65] and Es /∆ = 30, the fitted value is αc = 1.08. This indicates that the fit to (2.20) in Figure 12a gives an αc similar to those for other cutoff frequencies within numerical uncertainties. Such a finding was also reported in the previous NRG study16 where, even with three orders of magnitude change of ωc , the critical Kondo parameter remains essentially the same (Fig. 11 of reference 16.) This is physically incompatible with the prediction of the scaling limit αc → 1+ for ωc → ∞. IV.

CONCLUSIONS

In this paper we have used the multilayer improved relaxation approach51 to study the quantum phase transition in the spin-boson model. Such studies have been done with various approaches previously, the most relevant to our work is the NRG16 and the DMRG work.68 Different from our previous work51 employing equi-distance discretization of the phonon bath, here we use Wilson’s logarithmic discretization that is often found in NRG studies. This allows a faster (exponential) scaling of the bath frequency that can capture the infrared divergence more efficiently. Two quantities have been examined for quantification of the phase transition, the scaled energy splitting Es = ΛN (E1 − E0 ) and the scaled magnetic susceptibility Λ−N dhσz i/dǫ for the ground state with respect to the energy bias. The former is found to be more numerical appealing, although both provide similar information. Larger discrepancies were found, however, in how to utilize Es . A method that is often applied in NRG studies is to fit αc 15



Page 16 of 41

via (2.20) for α’s smaller than the critical value. However, we found such a procedure not reliable both numerically and physically. In particular, it misses the scaling limit αc → 1+ for ωc → ∞. To our knowledge there is no rigorous justification for using Eq. (2.20) to extract the critical coupling strength for the quantum phase transition. However, there is a speculation that Eq. (2.20) is valid in the asymptotic limit of α → αc , although we are not aware of any fitting of Eq. (2.20) in this limit numerically. To test this ansatz, we performed additional calculations in the range of α ∈ [0.91, 0.995] for several characteristic frequencies considered in this paper. Figure 13 shows the result for ωc /∆ = 5. The value of αc differs from that in Figure 10a, which is similar to what we have discussed in the previous section (that αc depends on the range of fitting α’s). More importantly, as we change the fitting range of α from [0.91, 0.995] in Figure 13a to [0.94, 1.03], the fitted αc also becomes larger and further away from the value of 1.08-1.1 obtained by examining the qualitative change of Es = ΛN (E1 − E0 ) before and after αc (Fig. 10b). Even more striking discrepancies are found in Fig. 14, where higher characteristic frequencies ωc are considered. The fitted αc values remain almost constant as ωc increases, and are much larger than the values extracted from examining the qualitative change of Es = ΛN (E1 − E0 ), presented in table 1. (Also see the NRG work, Fig. 11 of reference 16.) Furthermore, it is seen clearly from Fig. 14 that Eq. (2.20) is no longer a good description in the range α ∈ [0.91, 0.995]. While Eqs. (2.19)-(2.20) have been used extensively to extrapolate the critical coupling for the phase transition in the spin-boson model, we have found a significant discrepancy between this procedure and the examination of the qualitative change of Es = ΛN (E1 − E0 ) before and after αc . The values obtained from fitting (2.20) appear to be unreliable and do not possess the correct physical limit. Therefore, our findings suggest that more scrutiny is needed in such analysis. On the other hand, it was noted before70 that Eq.(2.20) breaks down when α is very close to αc (citation 38 in Ref. 16, although it was not considered as a problem.) In this case the fitting formula should be replaced by √

T ∗ ∝ econst/

16

αc −α

,


(4.1a)

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


or lnT ∗ = C + √

B , αc − α

(4.1b)

where C and B are constants. Figure 15 illustrates the fitting results using this expression. It is seen that the fitted αc values are in much better agreement with those extracted from examining the qualitative change of Es = ΛN (E1 − E0 ) as shown in table 1 (second and third columns.) Specifically, αc ’s for ωc /∆ = 5 and ωc /∆ = 10 agree with the lower bounds extracted from examining the qualitative change of Es = ΛN (E1 − E0 ), second column in table 1. The αc value in Figure 15c, for ωc /∆ = 20, is smaller than the lower bound in the corresponding second column of table 1. Our analysis shows that to correct this the range of α for fitting Eq. (4.1) needs to be brought even closer to αc . This is also true for ωc /∆ = 5 and ωc /∆ = 10. An approximate way to make corrections to the fit without performing additional calculations is by noticing that the least-square fit to Eq. (4.1) is unbiased to all α’s. On the logarithmic scale the plot shows a bigger deviation for α → αc where T ∗ becomes much smaller. This is undesirable because Eq. (4.1) is only valid in the α → αc limit. One may compensate this by multiplying a weighting function W (α) to bias the fit. As an example, we use the following weighting function W (α) = e(α−αmin )/δ ,

(4.2)

where αmin is the minimum α in the fitting data and δ is a strength parameter. Figure 16 shows such fittings with δ = 0.001, where the larger α’s (those closer to αc ) have higher weights. The fitted αc ’s for ωc /∆ = 5 and ωc /∆ = 10 now agree with the upper bound extracted from examining the qualitative change of Es = ΛN (E1 − E0 ), third column in table 1. The αc value in Figure 16c, for ωc /∆ = 20, is still a bit smaller than the lower bound in table 1. This suggests that more calculations are needed to bring the range of α even closer to αc , e.g., α ∈ [0.98, 0.9995]. Also, the multiplied weighting function is too sharp. Although its intention is to bias the fitting to α’s closer to αc , there may not be enough number of points in this range to warrant the accuracy of fitting. In summary, our finding reveals that using Eq. (4.1), the fitted critical coupling αc is in good agreement with the value extracted from examining the qualitative change of Es = ΛN (E1 − E0 ). Thus, one should use Eq. (4.1) instead of Eq. (2.20) to extract αc . A more robust procedure is to first estimate αc by examining the qualitative change of Es = 17



Page 18 of 41

ΛN (E1 − E0 ), and performs careful calculations for a range of α’s that are close to this estimated αc . These data are then fit to Eq. (4.1), including possibly a weighting function to bias the fit to α’s closer to αc . This procedure may be iterated several times to further improve the accuracy of fitting.

ACKNOWLEDGMENTS

We are grateful to an anonymous reviewer for insightful comments that lead to the detailed analysis presented in the Conclusions. HW acknowledges the support from the National Science Foundation CHE-1500285, and the 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. JS acknowledges the support from the National Natural Science Foundation of China (No. 21421003).

18


Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES

1

Marcus, R. A.; Sutin, N. Electron Transfers in Chemistry and Biology; Biochim. Biophys. Acta 1985, 811, 265–322.

2

Su´ arez, A.; Silbey, R. Properties of a Macroscopic System as a Thermal Bath; J. Chem. Phys. 1991, 95, 9115–9121.

3

Weiss, U.; Grabert, H.; Linkwitz, S. Influence of Friction and Temperature on Coherent Quantum Tunneling; J. Low Temp. Phys. 1987, 68, 213–244.

4

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.

5

Weiss, U. Quantum Dissipative Systems; World Scientific: Singapore, 2nd ed., 1999.

6

Harris, R. A.; Silbey, R. On the Stabilization of Optical Isomers through Tunneling Friction; J. Chem. Phys. 1983, 78, 7330–7333.

7

Garg, A.; Onuchic, J. N.; Ambegaokar, V. Effect of Friction on Electron Transfer in Biomolecules; J. Chem. Phys. 1985, 83, 4491–4503.

8

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.

9

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.

10

Lehle, H.; Ankerhold, J. Low Temperature Electron Transfer in Strongly Condensed Phase; J. Chem. Phys. 2004, 120, 1436–1449.

11

Egger, R.; Weiss, U. Quantum Monte Carlo Simulation of the Dynamics of the Spin-Boson Model; Z. Phys. B 1992, 89, 97–107.

12

Egger, R.; Mak, C. H. Low-Temperature Dynamical Simulation of Spin-Boson Systems; Phys. Rev. B 1994, 50, 15210–15220.

13

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.

14

Sim, E.; Makri, N. Path Integral Simulation of Charge Transfer Dynamics in Photosynthetic Reaction Centers; J. Phys. Chem. B 1997, 101, 5446–5458.

19



15

Page 20 of 41

M¨ uhlbacher, L.; Egger, R. Crossover from Nonadiabatic to Adiabatic Electron Transfer Reactions: Multilevel Blocking Monte Carlo Simulations; J. Chem. Phys. 2003, 118, 179–191.

16

Bulla, R.; Lee, H.-J.; Tong, N.-G.; Vojta, M. Numerical Renormalization Group for Quantum Impurities in a Bosonic Bath; Phys. Rev. B 2005, 71, 045122.

17

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.

18

Anders, F. B.; Bulla, R.; Vojta, M. Equilibrium and Nonequilibrium Dynamics of the SubOhmic Spin-Boson Model; Phys. Rev. Lett. 2007, 98, 210402.

19

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–114.

20

Stockburger, J. T.; Grabert, H. Exact c-Number Representation of Non-Markovian Quantum Dissipation; Phys. Rev. Lett. 2002, 88, 170407.

21

Shao, J. Decoupling Quantum Dissipation Interaction via Stochastic Fields; J. Chem. Phys. 2004, 120, 5053–5056.

22

Yan, Y.-A.; Shao, J. Stochastic Description of Quantum Brownian Dynamics; Front. Phys. 2016, 11, 110309.

23

Str¨ umpfer, J.; Schulten, K. Open Quantum Dynamics Calculations with the Hierarchy Equations of Motion on Parallel Computers; J. Chem. Theor. Comput. 2012, 8, 2808–2816.

24

Gelin, M.; Tanimura, Y.; Domcke, W. Simulation of Femtosecond Double-Slit Experiments for a Chromophore in a Dissipative Environment; J. Chem. Phys. 2013, 139, 214302.

25

Zheng, X.; Jin, J.; Welack, S.; Luo, M.; Yan, Y. Numerical Approach to Time-Dependent Quantum Transport and Dynamical Kondo Transition; J. Chem. Phys. 2009, 130, 164708.

26

Wang, H.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory; J. Chem. Phys. 2003, 119, 1289–1299.

27

Wang, H.; Thoss, M. Numerically Exact Quantum Dynamics for Indistinguishable Particles: The Multilayer Multiconfiguration Time-Dependent Hartree Theory in Second Quantization Representation; J. Chem. Phys. 2009, 131, 024114.

28

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.

29

Wang, H. Multilayer Multiconfiguration Time-Dependent Hartree Theory; J. Phys. Chem. A 2015, 119, 7951–7965.

20


Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


30

Frenkel, J. Wave Mechanics; Clarendon Press: Oxford, 1934.

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 Multiconfiguration Time-Dependent Hartree Approach; J. Chem. Phys. 2006, 124, 034114.

34

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.

35

Wang, H.; Thoss, M. Quantum Dynamical Simulation of Electron-Transfer Reactions in an Anharmonic Environment; J. Phys. Chem. A 2007, 111, 10369–10375.

36

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.

37

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.

38

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.

39

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.

40

Wang, H.; Thoss, M. From Coherent Motion to Localization: Dynamics of the Spin-Boson Model at Zero Temperature; New J. Phys. 2008, 10, 115005.

41

Egorova, D.; Gelin, M. F.; Thoss, M.; Wang, H.; Domcke, W. Effects of Intense Femtosecond Pumping on Ultrafast Electronic-Vibrational Dynamics in Molecular Systems with Relaxation; J. Chem. Phys. 2008, 129, 214303.

42

Velizhanin, K. A.; Wang, H.; Thoss, M. Heat Transport through Model Molecular Junctions:

21



Page 22 of 41

A Multilayer Multiconfiguration Time-Dependent Hartree Approach; Chem. Phys. Lett. 2008, 460, 325–330. 43

Velizhanin, K. A.; Wang, H. Dynamics of Electron Transfer Reactions in the Presence of ModeMixing: Comparison of a Generalized Master Equation Approach with the Numerically Exact Simulation; J. Chem. Phys. 2009, 131, 094109.

44

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.

45

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.

46

Albrecht, K. F.; Wang, H.; M¨ uhlbacher, L.; Thoss, M.; Komnik, A. Bistability Signatures in Nonequilibrium Charge Transport through Molecular Quantum Dots; Phys. Rev. B 2012, 86, 081412.

47

Wang, H.; Shao, S. Dynamics of a Two-Level System Coupled to a Bath of Spins; J. Chem. Phys. 2012, 137, 22A504.

48

Wang, H.; Thoss, M. Numerically Exact, Time-Dependent Study of Correlated Electron Transport in Model Molecular Junctions; J. Chem. Phys. 2013, 138, 134704.

49

Wang, H.; Thoss, M. Multilayer Multiconfiguration Time-Dependent Hartree Study of Vibrationally Coupled Electron Transport Using the Scattering-State Representation; J. Phys. Chem. A 2013, 117, 7431–7441.

50

Wilner, E. Y.; Wang, H.; Cohen, G.; Thoss, M.; Rabani, E. Bistability in a Nonequilibrium Quantum System with Electron-Phonon Interactions; Phys. Rev. B 2013, 88, 045137.

51

Wang, H. Iterative Calculation of Energy Eigenstates Employing the Multilayer Multiconfiguration Time-Dependent Hartree Theory; J. Phys. Chem. A 2014, 117, 9253–9261.

52

Wang, H.; Liu, X.; Liu, J. Accurate Calculation of Equilibrium Reduced Density Matrix for the System-Bath Model: A Multilayer Multiconfiguration Time-Dependent Hartree Approach and Its Comparison to a Multi-Electronic-State Path Integral Molecular Dynamics Approach; 31.

53

Meyer, H.-D.; Wang, H. On regularizing the MCTDH Equations of Motion; J. Chem. Phys. 2018, 148, 124105.

54

Wang, H.; Meyer, H.-D. On regularizing the ML-MCTDH Equations of Motion; J. Chem. Phys. 2018, 149, 044119.

22


Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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.; 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.

58

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.

59

Manthe, U.; Matzkies, F. Iterative Diagonalization within the Multi-Configurational TimeDependent Hartree Approach: Calculation of Vibrationally Excited States and Reaction Rates; Chem. Phys. Lett. 1996, 252, 71–76.

60

Matzkies, F.; Manthe, U. A Multi-Configurational Time-Dependent Hartree Approach to the Direct Calculation of Thermal Rate Constants; J. Chem. Phys. 1997, 106, 2646–2653.

61

Manthe, U. The State Averaged Multiconfigurational Time-dependent 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

Wang, H. Basis Set Approach to the Quantum Dissipative Dynamics: Application of the Multiconfiguration Time-Dependent Hartree Method to the Spin-Boson Problem; J. Chem. Phys. 2000, 113, 9948–9956.

64

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.

65

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.

66

Chakravarty, S. Quantum Fluctuations in the Tunnling between Superconductors; Phys. Rev. Lett. 1982, 49, 681–684.

67

Bray, A. J.; Moore, M. A. Influence of Dissipation on Quantum Coherence; Phys. Rev. Lett.

23



Page 24 of 41

1982, 49, 1545–1549. 68

Wong, H.; Chen, Z.-D. Density Matrix Renormalization Group Approach to the Spin-Boson Model; Phys. Rev. B 2008, 77, 174305.

69

Hur, K. L. Entanglement Entropy, Decoherence, and Quantum Phase Transitions of a Dissipative Two-Level System; Ann. Phys. 2008, 323, 2208–2240.

70

Kosterlitz, J. M. The Critical Properties of the Two-Dimensional XY Model; J. Phys. C: Solid State Phys. 1974, 7, 1046–1060.

24


Page 25 of 41

(a)

N

Λ (Ε1−Ε0)/∆

10 8 6

α = 0.3 α = 0.4 α = 0.5 α = 0.6 α = 0.7

4 2 100

50

N

150

200

150

200

(b) 1.1

Λ (Ε1−Ε0)/∆

α=1 α = 1.02 α = 1.035

1.0

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.9 0.8 50

100

N

FIG. 1: Dependence of the scaled energy splitting between the ground and the first excited states, ΛN (E1 − E0 )/∆, on the number of bath modes N . The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is Λ = 1.05.

25



(a)

-Λ d/d(ε/∆)

1 0.1

-N

0.01

α = 0.3 α = 0.4 α = 0.5 α = 0.6 α = 0.7

0.001 0.0001

100

50

N

150

200

(b) 1.2 1.1 1.0

-N

-Λ d/d(ε/∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 41

α=1 α = 1.02 α = 1.035

0.9 0.8

50

100

N

150

200

dhσz i FIG. 2: Dependence of the scaled derivative −Λ−N d(ǫ/∆) on the number of bath modes N for ǫ=0

the ground state. The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is

Λ = 1.05.

26


Page 27 of 41

(a) -1

10

-2

T*

10

-3

10

-4

10 1

2

1.5

2.5

1/(αc-α)

3

3.5

4

(b) -1

10 10

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-2

-3

10

10

-4

-5

10 1

2

3

1/(αc-α)

4

FIG. 3: The crossover scale T ∗ as a function of the Kondo parameter α ∈ [0.3, 0.8], fitted to Eq. (2.20). The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is Λ = 1.05. The values of cutoff scale and the critical coupling are (a) Es /∆ = 6, αc = 1.063; (b) Es /∆ = 7.5, αc = 1.049.

27



(a) -1

10

-2

T*

10

-3

10

10

-4

-5

10

2

3

4

1/(αc-α)

(b)

-1

10

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 41

10

-2

-3

10 1

1.5

2

2.5

1/(αc-α)

3

3.5

FIG. 4: The crossover scale T ∗ as a function of the Kondo parameter α ∈ [0.5, 0.85], fitted to Eq. (2.20). The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is Λ = 1.05. The values of cutoff scale and the critical coupling are: (a) Es /∆ = 5, αc = 1.080; (b) Es /∆ = 2.5, αc = 1.159.

28


Page 29 of 41

1.05 1.00

Λ (Ε1−Ε0)/∆

0.95

0.90

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.85

0.80 0.75

50

100

N

150

200

FIG. 5: Dependence of the scaled energy splitting ΛN (E1 − E0 )/∆ on the number of bath modes N for α ∈ [1.01, 1.035] with an increment of 0.0025. The result for α = 1.01 is the top curve and that for α = 1.035 is the bottom curve. The characteristic frequency of the bath is ωc /∆ = 10, and the scale factor is Λ = 1.05.

29



(a) -1

10

-2

T*

10

-3

10

-4

10 1

1.5

2

2.5

1/(αc-α)

3

3.5

4

(b) -1

10 10

-2

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 41

-3

10

10

-4

-5

10

2

3

1/(αc-α)

4

FIG. 6: The crossover scale T ∗ as a function of the Kondo parameter. The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is Λ = 1.01. The parameters appeared in Eq. (2.20) are: (a) α ∈ [0.3, 0.8], Es /∆ = 6, αc = 1.065; (b) α ∈ [0.5, 0.85], Es /∆ = 5, αc = 1.081.

30


Page 31 of 41

1.05 1.00

Λ (Ε1−Ε0)/∆

0.95

0.90

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.85

0.80 0.75

200

400

N

600

800

FIG. 7: Dependence of the scaled energy splitting ΛN (E1 − E0 )/∆ on the number of bath modes N for α ∈ [1.01, 1.035] with an increment of 0.0025. The result for α = 1.01 is the top curve and that for α = 1.035 is the bottom curve. The characteristic frequency of the bath is ωc /∆ = 10, and the scale factor is Λ = 1.01.

31



(a) -1

10

T*

10

-2

-3

10

10

-4

-5

10 1.5

2

3

2.5

3.5

1/(αc-α)

4

4.5

5

(b) 1.0

Λ (Ε1−Ε0)/∆

0.9

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 41

0.8 0.7 4

8

N

12

16

20

FIG. 8: Two means of extrapolating αc for the characteristic frequency of the bath is ωc /∆ = 10. The scaling factor is Λ = 2 for the logarithmic discretization. (a) The crossover scale T ∗ as a function of the Kondo parameter. The parameters appeared in Eq. (2.20) are: α ∈ [0.5, 0.85], Es /∆ = 5, αc = 1.067. (b) Dependence of the scaled energy splitting ΛN (E1 − E0 )/∆ on the number of bath modes N for α ∈ [1.01, 1.035] with an increment of 0.0025. The result for α = 1.01 is the top curve and that for α = 1.035 is the bottom curve. The estimated αc is 1.0275-1.03.

32


Page 33 of 41

1.0

-N

-Λ /(ε/∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.9 0.8

−4

ε/∆ = 10 −5 ε/∆ = 10 −7 ε/∆ = 10

0.7 0.6 0.5

50

100

N

150

200

FIG. 9: Dependence of the scaled numerical derivative −Λ−N hσz iǫ /(ǫ/∆) on the number of bath modes N for the ground state. The characteristic frequency of the bath is ωc /∆ = 10, and the scaling factor is Λ = 1.05.

33



(a)

-1

T*

10

10

-2

-3

10 1

2

1.5

3

3.5

800

1000

2.5

1/(αc-α)

(b) 1.8 1.6 1.4

Λ (Ε1−Ε0)/∆

1.2 1.0

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 41

0.8

0.6

0.4 0.2

200

400

600

N

FIG. 10: Two means of extrapolating αc for the characteristic frequency of the bath is ωc /∆ = 5. The scaling factor is Λ = 1.01. (a) The crossover scale T ∗ as a function of the Kondo parameter. The parameters appeared in Eq. (2.20) are: α ∈ [0.3, 0.8], Es /∆ = 3, αc = 1.085. (b) Dependence of the scaled energy splitting ΛN (E1 − E0 )/∆ on the number of bath modes N for α ∈ [1, 1.14] with an increment of 0.02. The result for α = 1 is the top curve and that for α = 1.14 is the bottom curve. The estimated αc is 1.08-1.10.

34


Page 35 of 41

(a) -1

10

T*

10

-2

-3

10

10

-4

-5

10 1

1.5

2

2.5

1/(αc-α)

3

3.5

4

(b) 1.05

Λ (Ε1−Ε0)/∆

1.00

0.95

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.90 200

400

600

N

800

1000


35



(a) -1

10

-2

T*

10

-3

10

10

-4

-5

10

1.5

2

1/(αc-α)

2.5

(b) 1.001 1.000

Λ (Ε1−Ε0)/∆

0.999

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 41

0.998 0.997

0.996

200

400

N

600

800


36


Page 37 of 41

(a) 10

-2

-3

T*

10

10

-4

-5

10 4.5

5

5.5

6

1/(αc-α)

6.5

7

7.5

(b) -1

10 10

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-2

-3

10

10

-4

-5

10 5

6

7

1/(αc-α)

8

9

FIG. 13: The crossover scale T ∗ as a function of the Kondo parameter α, fitted to Eq. (2.20). The characteristic frequency of the bath is ωc /∆ = 5, and the scaling factor is Λ = 1.01. The range of coupling, value of cutoff scale and critical coupling are: (a) α ∈ [0.91, 0.995], Es /∆ = 2, αc = 1.130; (b) α ∈ [0.94, 1.03], Es /∆ = 1.2, αc = 1.136.

37



(a) -1

10

-2

T*

10

-3

10

10

-4

-5

10

6

5

1/(αc-α)

7

8

(b) 10

-2

T*

10

0

10

-4

5.5

6

6.5

7

1/(αc-α)

7.5

8

8.5

(c) 0

10

10

-2

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 41

10

-4

7.5

8

1/(αc-α)

8.5

9

FIG. 14: The crossover scale T ∗ as a function of the Kondo parameter α ∈ [0.91, 0.995], fitted to Eq. (2.20). The scaling factor is Λ = 1.01. The values of frequency, cutoff scale, and the critical coupling are: (a) ωc /∆ = 10, Es /∆ = 1.23, αc = 1.117; (b) ωc /∆ = 20, Es /∆ = 1.09, αc = 1.117; (c) ωc /∆ = 100, Es /∆ = 1.054, αc = 1.105.

38


Page 39 of 41

(a)

-2

T*

10

10

-4

3.0

2.5

3.5

1/(αc-α)

1/2

4.0

4.5

(b)

T*

10

-2

-3

10

10

-4

-5

10 3.5

4.0

4.5

5.0

1/(αc-α)

5.5

1/2

6.0

6.5

(c)

10

-2

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10

-4

4

6

8

1/(αc-α)

10

1/2

12

FIG. 15: The crossover scale T ∗ as a function of the Kondo parameter α fitted to Eq. (4.1). The scaling factor is Λ = 1.01. The values of the fitting range, frequency, cutoff scale, and critical coupling are: (a) α ∈ [0.94, 1.03], ωc /∆ = 5, Es /∆ = 1.2, αc = 1.08; (b) α ∈ [0.95, 0.995], ωc /∆ = 10, Es /∆ = 1.23, αc = 1.02; (c) α ∈ [0.96, 0.995], ωc /∆ = 20, Es /∆ = 1.0925, αc = 1.002.

39



(a)

-2

T*

10

10

-4

3.0

2.5

4.0

3.5

1/(αc-α)

1/2

(b) -1

10

-2

T*

10

-3

10

10

-4

-5

10

4.0

3.5

4.5

1/(αc-α)

5.0

1/2

(c) 10

10

0

-2

T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 41

10

-4

4

5

6

7

8

1/(αc-α)

1/2

9

10

11

FIG. 16: Same as Fig. 15 but multiplying a weighting function in the form of Eq. (4.2) with δ = 0.001. The fitted critical coupling are: (a) αc = 1.092, (b) αc = 1.034; (c) αc = 1.0044.

40


Page 41 of 41

Table of Content Graphic

10

Frequency of the mode

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 0.1 0.01 0.001 0.0001 0

200

400

600

800

Number of discrete modes

41


1000

Quantum Phase Transition in the Spin-Boson Model: A Multilayer

Recommend Documents