From coherent motion to localization: dynamics of the spin-boson model at zero temperature

The dynamics of the spin-boson model at zero temperature is investigated using the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method. This method allows a numerically exact description of the dynamics of the spin-boson model in a broad range of coupling strengths. The results show the transition of the dynamics from weakly damped coherent motion to incoherent decay and finally to localization upon increase of the system–bath coupling strength. A detailed analysis of the incoherent decay for stronger coupling reveals that multiple timescales are involved in the dynamics. Furthermore, the applicability of the scaling limit for large characteristic frequencies of the bath as well as the validity of the non-interacting blip approximation are studied in some detail.


Introduction
The accurate description of quantum effects in dissipative systems is a central task in condensed matter physics. Due to the complexity of condensed phase systems, a variety of different models have been proposed, which focus on the important physical aspects of the problem under consideration and are amenable to a theoretical treatment. The spin-boson model [1]- [3] is a classical example in this regard for describing a variety of different processes including electron transfer [4], hydrogen tunneling [5], macroscopic quantum coherence [6], and many other phenomena [2] in the condensed phase. The model describes two states that are linearly coupled to a bath of harmonic oscillators. In mass-weighted coordinates, the Hamiltonian reads where σ x and σ z are Pauli matrices σ x = |φ 1 φ 2 | + |φ 2 φ 1 |, (1.2a) In the context of electron transfer theory [4], which we will mainly refer to in this work, the two states correspond to diabatic electronic states describing the donor (|φ 1 ) and the acceptor (|φ 2 ) state of the reaction, respectively. Depending on the specific application, the bath of harmonic oscillators models the solvent, the phonons of a solid, or other condensed phase environments. The coupling of the electronic states to the bath and all properties of the bath that influence the dynamics of the two-state subsystem are characterized by the spectral density [1,2] J B (ω) = π 2 j c 2 j 3 to describe the quantum dynamics of this model in appropriate physical regimes, including the non-interacting blip approximation (NIBA) [1], the Bloch-Redfield equation [1], [7]- [10], the master-equation/Smoluchowski-equation treatment of the reaction coordinate [11]- [16], mixed quantum classical and semiclassical methods [17]- [21], and other approaches [2,3].
On the other hand, numerical path integral calculations [22]- [26] based on Feynman-Vernon influence functional [27,28] have been performed to study the dynamics of the spin-boson model for weak to moderate system-bath coupling. Recently, we have reported a rather comprehensive study of the spin-boson model [29] employing the numerically exact selfconsistent hybrid approach [30]. In this study, we have also analyzed the validity of various approximate treatments in the different physical regimes. As pointed out in [29], there were at least two regimes where none of the tested approximate methods gave reliable results: (i) the regime with an intermediate electronic coupling (between adiabatic and nonadiabatic) and strong system-bath coupling and, (ii) the moderately adiabatic regime with strong system-bath coupling at low temperature.
In this paper, we present a study of the dynamics of the spin-boson model in a different physical regime-the zero temperature limit of the phonon bath. In contrast to our previous work on electron transfer reactions in chemical systems, which focuses on the dynamics of the spin-boson model at finite temperatures, the zero-temperature limit is of interest from the more fundamental physical point of view. In particular, the zero-temperature regime is characterized by strong quantum effects, a long memory time of the bath, and the transition from coherent damped motion to localization upon increase of the system-bath coupling strength. It is also known to be a challenging problem for numerical as well as analytical studies. Recently, there has been an impressive development of numerical simulation methods to treat the dynamics of this regime, in particular, the numerical renormalization group theory [31,32] and the mixed stochastic-deterministic scheme [33]. Our study in this paper is based on the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method [34]. This method allows us to obtain numerically exact results for the dynamics of the spin-boson model at zero temperature in a broad range of coupling parameters. Due to the rigorous treatment, the ML-MCTDH method also provides benchmark results to test approximate approaches.
The remainder of the paper is organized as follows. In section 2, we outline the basic idea of the ML-MCTDH method as well as some details specific to the application to the spin-boson system at zero temperature. In section 3, we present results obtained from our numerically exact simulations. In particular, we discuss the transition from weakly damped coherent motion to incoherent decay and finally to localization when the system-bath coupling strength is increased. Furthermore, we discuss the validity of the approximate NIBA approach and the applicability of the scaling limit. Section 4 summarizes our findings.

Observables of interest, initial state and discretization of the bath
Various observables are of interest to study the dynamics of the spin-boson system, e.g. the population or the coherence of the two states, the energy of the bath oscillators, or linear and nonlinear optical spectra. In this paper, we will focus specifically on the time-dependent 4 population difference of the two level system where we have used atomic units in whichh = 1. Thereby, we consider a factorized initial state specified by the (donor) state |φ 1 for the two level system and the Boltzmann operator e −β H B for the bath. For a theoretical study of spin-boson type models with a correlated initial state and time-dependent laser fields using the ML-MCTDH method see [35].
For the specific case of zero temperature, considered in this work, the initial state becomes a pure state-the ground state |χ 0 of H B . Accordingly, equation (2.1) is rewritten as In most of the results presented below, H B is chosen as describing a nonequilibrium initial preparation that can be achieved, e.g. by photoexcitation from a lower-lying electronic state with equilibrium geometry between the donor |φ 1 and the acceptor state |φ 2 . Other initial conditions considered below are where H L B and H R B correspond to a bath in equilibrium with the charge distribution of the donor and the acceptor state, respectively.
As discussed in the introduction, the bath and its coupling to the two-level system is fully characterized by the spectral density J (ω) in equation (1.3). In the applications considered below, we will consider a spectral density of Ohmic form with exponential cutoff 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 spin-boson model to the Kondo problem has been discussed previously [1,2]. In the context of electron transfer theory, the system-bath coupling is often specified by the reorganization energy λ, which is related to the Kondo parameter and the characteristic frequency of the bath via λ = 2παω c .
Other, more complex forms of the spectral density as well as an anharmonic bath have been studied previously at finite temperatures [36]- [38]. The zero-temperature dynamics of these models will be reported in future publications. The spectral density (2.4) describes a condensed phase environment with a formally continuous distribution of bath modes. Within the timescale of interest, a finite number of bath modes can be used to represent the condensed phase environment [29,30,34]. Similar discretization methods have been used in many other numerical procedures, e.g. Gaussian quadrature methods, discrete Fourier transforms and the finite element method. The number 5 of bath modes required to obtain a proper representation of the condensed phase environment depends on the timescale of interest and the specific model parameters. To ensure convergence to the condensed phase limit, we have employed a few hundred to a few thousand modes for the various cases discussed in this paper.
The continuous bath spectral density of equation (2.4) can be discretized to the form of equation (1.3) via the relation [29,30,34] where ρ(ω) is a density of frequencies satisfying In this work, ρ(ω) is chosen as where N b is the number of discrete bath modes in the simulation.
To improve the efficiency of the calculations, some of the high-frequency modes can be removed from the dynamical simulation using a Born-Oppenheimer type approximation [1], resulting in a modified electronic coupling Here, ω q defines the boundary frequency above which the modes are not treated explicitly in the dynamical simulation but are included via the effective coupling parameter eff . This procedure has also been used in the framework of the self-consistent hybrid approach [29,30,39]. As for the number of modes used in the discretization procedure discussed above, the boundary frequency ω q is a numerical convergence parameter. Consequently, the actual number of modes which can be treated via equation (2.7), without deteriorating the accuracy of the calculation, depends on the specific parameters. Our results show that in situations with weak system-bath coupling (i.e. small α) one can generally remove a large percentage of the high-frequency modes from the dynamical calculation and treat them via equation (2.7). As α becomes larger, converged results require that fewer modes are treated this way. For the strongest coupling considered in this paper (α = 1.5), no modes can be treated via equation (2.7).

Multilayer multiconfiguration time-dependent Hartree theory
To obtain the time-dependent population difference P(t) in equation (2.2), the major task is the time evolution of the wave function for a system including many degrees of freedom. This is achieved by employing the recently developed ML-MCTDH theory [34]. The method as well as its applications to reactions in the condensed phase have been described in detail previously [34,35,38], [40]- [42]. Here, we only briefly introduce the general idea and give some details specific to the application in this work. In addition, we outline explicitly the ML-MCTDH theory for more than two dynamical layers, which is crucial for the applications considered here.

6
The ML-MCTDH method [34] is a rigorous variational approach to study quantum dynamics in systems with many degrees of freedom. It extends the original MCTDH method [43]- [46] for applications to significantly larger systems. In the original (singlelayer) MCTDH method, the overall wave function is expanded in terms of time-dependent configurations Each configuration is a Hartree product (for systems with distinguishable particles) of the 'single-particle' functions (SPFs) |ϕ k j k (t) , where M denotes the total number of single particle (SP) degrees of freedom and k is the index of a particular SP group. In practice, one SP group usually contains several (Cartesian) degrees of freedom. In the MCTDH method, the SPFs are represented by the full configuration-interaction (CI) expansion of the time-independent basis functions, This limits the application of the MCTDH method to a few tens of degrees of freedom [43]- [46]. In the ML-MCTDH theory, the full CI construction of the SPFs in equation (2.9) is replaced by a time-dependent multiconfigurational expansion (2.10) Here, Q(κ) denotes the number of level two (L2) SP groups in the κth level one (L1) SP and |v (κ,q) i q (t) is the L2-SPF for the qth L2-SP group. Both are contained in the κth L1-SP group. The expansion of the overall wave function can be formally written as can be expressed by a further multiconfigurational expansion. As a result, the overall wave function | (t) can be expanded recursively to many layers in the ML-MCTDH framework. Following the Dirac-Frenkel variational principle [47], the equations of motion are obtained from variation of the wave function | (t) with respect to the expansion coefficients of each layer (t) are reduced density matrices for the first, second and third layers, respectively, and Ĥ (κ) (t), Ĥ (κ,q) L2 (t) and Ĥ (κ,q,γ )

L3
(t) are mean-field operators for the corresponding layers. These operators can be recursively evaluated by means of the single hole functions L3;m,s (t) for the first, second, third, and further layers.
The single hole function for the first layer is defined by projecting a particular SPF out of the overall wave function Hole functions for other layers are defined in a similar way by projecting a particular SPF for this layer out of the SPF for the upper layer. For example, the single hole function |g (κ,q) L2;n,r (t) for the second layer can be constructed from the relation (t), . . . , are the SP-space projection operators for different layers [34].
The inclusion of several dynamically optimized layers in the ML-MCTDH method provides more flexibility in the variational functional, which significantly advances the capabilities of performing wave packet propagations in complex system. This has been demonstrated by several applications to quantum dynamics in the condensed phase including many degrees of freedom [34,35,37,40,41]. In this work, we employ an implementation of the ML-MCTDH theory with up to four dynamical layers plus one static layer. The number of degrees of freedom in the simulation ranges from a few hundred to a few thousand. Unless specified otherwise, only converged quantum results are shown. Readers interested in the numerical details of the method and the convergence procedure are referred to our earlier published work [34,35,37,40,41].

Results and discussion
The methodology outlined above is applied to study the dynamics of the unbiased ( = 0) spin-boson model with Ohmic spectral density at zero temperature. In particular, we shall consider the transition from weakly damped coherent motion to incoherent decay and finally to localization when the system-bath coupling strength is increased. Using the numerically exact results of the ML-MCTDH method, we will also discuss the validity of the NIBA approach.

Weak to moderate coupling regime: from coherent motion to incoherent decay
We first consider the weak to moderate coupling regime (α 0.5) with a nonequilibrium initial condition defined by H B = H 0 B in equation (2.3a). Figure 1 shows the time-dependent population difference, equation (2.2), for different characteristic frequencies of the bath: (a) ω c / = 10, (b) ω c / = 20, (c) ω c / = 40, and different Kondo parameters, α = 0.05-0.5. It is instructive to analyze the influence of the system-bath coupling and the characteristic frequency of the bath by varying only one of the two parameters at a time. Considering results for fixed ω c , it is seen from figure 1 that an increase of the Kondo parameter α, which determines the system-bath coupling strength, introduces (electronic) decoherence effects in the dynamics of 9 P(t). Specifically, when α increases from 0.05 to 0.5, the characteristics of the population dynamics changes from weakly damped coherent motion to incoherent decay. On the other hand, considering results for the same value of α, increasing ω c also influences the dynamics of P(t) significantly, both in the coherent and incoherent regimes. In the coherent regime, the period of the coherent oscillation becomes longer and coherence is less pronounced when ω c increases. In the incoherent regime, the timescale of the incoherent decay becomes slower as ω c increases.
The above findings are consistent with electron transfer theory [4] in which the reorganization energy is a measure for the overall electronic-nuclear coupling For the Ohmic spectral density, λ = 2παω c is proportional to both the Kondo parameter α and the characteristic frequency ω c . Thus, the increase in either quantity will reduce electronic coherence or increase the timescale for the incoherent decay.
It has been argued that in the so-called 'scaling limit' of a fast bath (i.e. ω c / 1), the characteristic frequency of the environment, ω c , influences the population dynamics only through a scaling of the time by the following renormalized coupling parameter [2] .
(3.2) Figure 2 illustrates the dependence of P(t) on the scaled time t r for three Kondo parameters, α = 0.05, 0.2 and 0.5. It is seen that as ω c increases, P(t) approaches the scaling limit as discussed in previous work [2]. For the Ohmic spectral density (with exponential cutoff) considered in this paper, the scaling limit is approximately reached when ω c / ≈ 20-40, corresponding to a very fast timescale of the bath. In accordance with equation (3.2), the results indicate that in this limit the oscillation periods are proportional to ω α/(1−α) c . Similarly, the timescale of the incoherent relaxation agrees well in the scaled time t r (cf figure 2(c)) indicating that the original decay becomes slower with increasing ω c .
However, figure 2 also shows that for moderately large values of the characteristic frequency of the bath, ω c / ≈ 5-10, P(t) deviates appreciably from the results for larger ω c . Thus, in this parameter regime, the influence of ω c on the dynamics of P(t) is more complicated than a simple scaling of the time. In particular, the results demonstrate clearly that an increase of ω c results in a quenching of coherence, which is not captured by the scaling law (3.2). On the other hand, the period of the oscillations in the weakly damped regime, figures 2(a) and (b), follows the scaling law.
It is noted that in this paper, we focus on the dynamics of the spin-boson model for a relatively fast bath (ω c / > 1). The influence of ω c in a broader range has been discussed previously [29] for the spin-boson model at finite temperatures. An analysis of the zero temperature dynamics will be presented in the future.
We now discuss the performance of approximate NIBA approach in the weak to moderate coupling regime (α 0.5). between the NIBA results and those from the numerically exact simulation are better for larger ω c as compared with the corresponding results in figure 3. In accordance with previous results [2,33,36], these findings demonstrate that NIBA is a good approximation in the weak coupling (α 1), nonadiabatic (ω c / 1) regime. Although, in general, NIBA becomes less accurate when α increases, there is a special case of Kondo parameter α = 0.5, where NIBA gives the 'exact' result in the scaling limit ω c → ∞ [1,2]. In this limit and for zero temperature, a simple expression for P(t) can be obtained (cf for example, equation (5.23) in [1], note there is a factor of two difference in the convention of between [1] and the current paper) Therefore, α = 0.5 is sometimes loosely referred to as the 'exact' regime for NIBA. In agreement with the discussion above, figure 5 shows that the applicability of the scaling limit depends on how large ω c is compared to the nonadiabatic coupling . For ω c / = 10, the result of NIBA shows deviations from the numerically exact result. The agreement becomes better as ω c / increases, similar to the comparisons discussed above for other values of α.
When ω c / = 40, NIBA becomes almost exact, indicating that the scaling limit is reached for the Ohmic spectral density. The results shown above have been obtained for a nonequilibrium preparation of the initial state of the bath described by equation (2.3a). Figure 6 illustrates the dependence of the population dynamics, P(t), on different initial conditions defined in equation (2.3). In the weak coupling regime, α = 0.05, virtually no influence of the bath initial condition on P(t) can be found (cf figure 6(a)). Small differences, mainly within the first oscillation period, are observed when α is increased to 0.2 ( figure 6(b)). For a stronger coupling strength, α = 0.5, the time-dependent population eventually shows appreciable dependence on the bath initial condition. Although NIBA can qualitatively predict this trend (data not shown), the results deviate from the numerically exact simulation when the coupling strength becomes large, as shown in figures 3-5.

Moderate to strong coupling regime: characteristics of the incoherent decay
We next consider the dynamics of the spin-boson model in the moderate to strong coupling regime, 0.5 < α < 1. More precisely, the lower boundary of this parameter regime is defined by the coherent to incoherent transition (α ∼ 0.5 in the scaling limit) and the upper boundary by the transition to localization (α loc 1, see the next subsection). In this regime, the dynamics of P(t) is fully incoherent. The numerically exact results depicted in figure 7 show that the decay of P(t) becomes slower for increasing α. Furthermore, for fixed Kondo parameter α, P(t) decays slower with increasing characteristic frequency of the bath, ω c . This qualitative trend has also been found in previous studies [2,23,32,33]. NIBA also shows this trend but predicts a qualitatively incorrect purely algebraic decay [2]. The analysis in the following subsection reveals that the transition to localization occurs at α loc 1.2 for the parameters used in figure 7(a), whereas for figure 7(c) α loc is much closer to one. As a consequence, the value α = 0.9 is much further away from the transition to localization in figure 7(a) as compared to figure 7(c). Similar conclusions are true for other values of α.
We now turn to a more quantitative analysis of the population dynamics in the parameter regime 0.5 < α < 1. Motivated by the study in [33], we first analyzed the dependence of ln[P(t)] versus t to identify a (possible) long-time exponential decay. This analysis has been done for all simulation results in this parameter regime. As a representative example, figure 8 shows ln[P(t)] ∼ t plots for two different coupling strengths and several different characteristic frequencies. At first glance, the results shown in figure 8 seem to indicate that the long-time characteristics is indeed exponential decay, i.e. ln[P(t)] depends linearly on t. However, a quantitative analysis employing a fit of the expression ln[P(t)] = −kt after a 'sufficiently long' cutoff time t 0 reveals that this is not the case. Table 1 summarizes the results for the fitted exponential rate constant k. The results indicate that the value of k depends sensitively on the cutoff time t 0 . Specifically, the values obtained for k decrease with increasing t 0 . One may argue that in the limit of t 0 → ∞ a converged rate constant k can be obtained. However, such a rate constant would be of little relevance to characterize the dynamics of P(t), especially for 0.5 < α < 0.75, because (as shown in figures 7 and 8) a significant decay of P(t) has already occurred within the time t 0 = 10. Furthermore, the cutoff time t 0 should represent the period within which P(t) undergoes partially coherent, transient dynamics, before reaching incoherent decay. This transient timescale is expected to be long for small Kondo parameters α and to decrease with increasing α. When α is sufficiently large (i.e. larger than 0.5), t 0 is expected to be short and of the order of t 0 ∼ 1. Instead, our single exponential fit requires t 0 to be large even for very large α and, furthermore, t 0 increases versus α. Thus, a more flexible ansatz has to be used to characterize the dynamics of P(t). We have found that a multi-exponential fit with two or three terms provides an excellent description for P(t) over the entire range of simulation (t < 80). In this procedure P(t) is fitted, after a short Thereby, since the fitting starts at t 0 , i a i = 1. This is correct because P(t) always starts as a concave, non-exponential curve, and changes to incoherent decay only after the transient time t 0 . The fact that a relatively short timescale of t 0 = 1 is sufficient to account for the transient behavior shows that the physics of the problem is described correctly. For 0.5 < α < 0.75, the fitting is performed for t 40 (or shorter time if P(t) is less than 0.01). For α 0.75, the fitting range is t 80. The fact that a multi-exponential fit of the form of equation (3.4) is required implies that multiple timescales are involved in the dynamics of P(t). We note that the general form of equation (3.4) has been predicted by conformal field theory (CFT) [48]. However, the thus obtained decay rates k i are quite different. In contrast to the results of CFT, the rate constants obtained by the multi-exponential fit, equation (3.4), are not integer multiples of a single rate. On the other hand, analysis of the (approximate) dynamics at asymptotic times suggested that P(t) might be described by an exponential or a combination of exponential and algebraic decay [2,48,49]. This asymptotic behavior, however, may only be realized at a sufficiently long time t r 1, where r is defined in equation (3.2). As discussed below with respect to the scaling limit (also see figure 10 for relevant timescales), such asymptotic times are orders of magnitude longer than the timescales that are relevant for the relaxation dynamics of P(t) and will not be considered here.   Tables 2-4 show results of fitting P(t) to a bi-exponential decay form, i = 1, 2 in equation (3.4). Over the entire range of 0.5 α 1, this bi-exponential fit gives quite an accurate description of the dynamics of P(t). This is illustrated in figure 9. Thus, in general, there are two timescales that dominate the decay dynamics of P(t): a fast exponential decay (characterized by k 1 ) that is responsible for the short time characteristics (after the transient time t 0 ), and a slower decay (characterized by k 2 ) that describes the long-time incoherent dynamics of P(t). For 0.5 < α < 0.75, both k 1 and k 2 decrease with increasing Kondo parameter α. For even larger values of α the bi-exponential is less accurate (see below) and thus the trend is blurred. Comparing tables 2-4 with table 1, it is seen that k 2 approximately corresponds to the decay constant obtained from the fit to the single exponential form in the asymptotic region.
For coupling strengths α larger than 0.75, the bi-exponential fit to the time-dependent population becomes less accurate. Table 5 displays results obtained using a tri-exponential fit to P(t) in the range of 0.75 α 1. Excellent agreement is observed between the fit and the simulation results, as shown in figure 9. This result indicates that more and more timescales are involved in the dynamics if the system-bath coupling strength α approaches unity, which could also imply that the decay becomes non-exponential. Figure 9 also shows that the single Table 3. The bi-exponential fit to P(t), equation (3.4), for ω c / = 20. The fitting starts at t 0 = 1, and the constants k i are reported as k i / . To provide more physical insight into the dynamics of P(t) for the moderate to strong coupling regime (0.5 < α < 1), we have analyzed the population dynamics P(t) as a function of the scaled time t r , where r is given in equation (3.2), similar to the study of the scaling limit for the coherent regime, α 0.5, in figure 2. In contrast to situations where α 0.5, the exponent in equation (3.2), α/(1 − α), is larger than unity in the strong coupling regime. As a result, r changes more rapidly with respect to variations of α or ω c . Therefore, the realization of the scaling limit requires stricter physical and mathematical criteria. For example, equation (3.2) predicts that for α = 1 ( r = 0) the dynamics of P(t) is characterized by localization in the scaling limit ω c → ∞. However, with a finite value of ω c , the population still undergoes a slow decay (see below). In this case equation (3.2) cannot be used to predict the change of P(t) versus ω c . Similarly, when the value α is close to unity, the exponent α/(1 − α) is very large, which requires a very long timescale to examine the scaling limit. Thus, for the analysis of the scaling limit, we focus on the range of 0.5 < α 0.75, which is tractable by time-dependent studies. Figure 10 illustrates (similar to figure 2) the dependence of P(t) on the scaled time t r for three Kondo parameters, α = 0.55, 0.65 and 0.75. Results for other values of α in this range show similar qualitative behavior (data not shown). When α is close to 0.5, i.e. the exponent in equation (3.2) is close to unity, only the results with smaller characteristic frequency ω c / = 10 show deviations from results for larger ω c (cf figure 10(a)). Thus, in this regime, one may conclude that the scaling limit is reached when ω c / ≈ 20 − 40, similar to the results in figure 2 for α < 0.5. However, as the Kondo parameter α increases, the realization of the scaling limit requires a larger characteristic frequency. For example, when α = 0.65-0.75, figures 10(b) and (c) show that the scaling limit is approximately reached with ω c / ≈ 40-60. This finding is consistent with previous analytic work [2]. Recent numerical simulations [33], on the other hand, seem to indicate that the scaled time t s is independent of the Kondo parameter, s = ( /ω c ). We note, however, that the form of the spectral density in [33] is different from the present work and analytic theories [1,2]. In addition, Zhou and Shao [33] investigated the long-time scaling limit of the rate constant with respect to s , whereas the current work focuses on the characteristic of P(t) in the entire range of the simulation time. The present work shows that P(t) obeys the scaling behavior in the appropriate limit even if it cannot be described by the single exponential decay.

Strong coupling regime: from incoherent decay to localization
With further increase of the Kondo parameter α, the decay of P(t) becomes slower and eventually reaches localization. This transition from incoherent decay to localization is illustrated in figure 11(a) for ω c / = 10. The specific coupling strength α loc above which localization occurs depends on the characteristic frequency of the bath, ω c . For ω c / = 10, the transition to localization takes place at α loc = 1.2-1.3. As ω c increases, α loc gets closer to unity.
In the limit ω c / → ∞, the transition to localization occurs at α loc = 1. This transition has been predicted by various analytical theories [1,2,50,51], rigorous mathematical arguments [52], the flow-equation approach [53] as well as NRG calculations [31,32]. The extent of localization, as measured by the deviation of P(∞) from its thermal equilibrium value 1/2, also depends on ω c . Results of simulations for different values of ω c , depicted in figure 11(b), indicate that in the limit ω c → ∞ maximal localization is obtained, i.e. lim ω c →∞ P(∞) = 1.

Computational effort
To conclude this section, we briefly comment on the numerical effort required to obtain the above simulation results. For weak to moderate coupling strength (0 < α < 0.5), a few hundred bath modes provide numerically exact results for the dissipative dynamics within t = 80. In this regime, converged ML-MCTDH results take between 10 min and a few hours on a standard Pentium 4 (3 GHz) computer. The stronger coupling regime (0.5 < α < 1) is more challenging, where 500-1000 bath modes are required to achieve convergence. In this regime, each converged result takes about 8-16 h on a Pentium 4 computer. For even stronger coupling strength (α > 1), the simulation becomes more demanding: 1000-3000 modes are required to represent the condensed phase environment properly. The basis function for each mode ranges from three (for high-frequency modes) to a few hundreds (for low frequency modes). The configuration space for each layer [34] is typically a few hundred thousands. This large variational space is caused by the strong quantum effects in this regime and the fact that we require the final results be converged to within 3-4 significant figures (cf figure 11). In the strong coupling regime, each converged simulation took a few days to a week. On the other hand, we are unaware of other methods that can achieve such a high accuracy in this interesting though difficult physical regime. The computational costs described above refer to converged quantum dynamical calculations up to a time of t = 80. To ensure that convergence is achieved, for each physical parameter a series of careful convergence tests are performed with respect to all the variational parameters such as the number of bath modes, primitive basis functions, SPFs in each layer, etc. Based on these tests, we estimate that the results presented in this paper are converged to within a few per cent relative error within the timescale of t = 80. We note that for even longer time more bath modes are necessary to converge the dynamics. A few such tests have been carried out and revealed that none of the qualitative conclusions made above would be altered, suggesting that a sufficiently long timescale of the dynamics of the spin-boson model has been explored in this work.

Concluding remarks
In this paper, we have studied the dynamics of the spin-boson model at zero temperature. This is a particularly interesting and challenging parameter regime characterized by strong quantum effects and a long memory time of the bath. To simulate the dynamics of the spin-boson model in this parameter regime, we have used the ML-MCTDH method, which provides numerically exact results for a broad range of coupling strengths.
The results show the transition of the dynamics from weakly damped coherent motion to incoherent decay and finally to localization upon increase of the system-bath coupling strength. A detailed analysis of the decay characteristics for stronger system-bath coupling demonstrates, contrary to previous studies, that for relevant times the dynamics cannot be described by a single exponential rate constant but involves more timescales.
We have also analyzed the validity of the approximate NIBA approach and the applicability of the scaling limit, where the characteristic frequency of the bath influences the dynamics of the two state system only via a scaled time variable. The simulations show that this limit is realized only for rather large characteristic frequencies of the bath, the actual value depends on the system-bath coupling strength.
In this work, we have focused on the unbiased spin-boson model with Ohmic spectral density. The study of other, more complex spectral densities [36,37] as well as the influence of an energy bias will be the subject of future work.