Open access peer-reviewed chapter

Unitary Approaches to Dissipative Quantum Dynamics

Written By

Matteo Bonfanti and Rocco Martinazzo

Submitted: 30 October 2015 Reviewed: 24 February 2016 Published: 24 August 2016

DOI: 10.5772/62686

From the Edited Volume

Research Advances in Quantum Dynamics

Edited by Paul Bracken

Chapter metrics overview

1,741 Chapter Downloads

View Full Metrics


We describe in detail a general system–bath strategy for investigating the quantum behavior of small systems interacting with complex environments. In this approach, a simplified heat bath is used as a surrogate for realistic environments, and explicit, unitary quantum simulations of the “universe” (the system plus the bath) are performed by means of high-dimensional wave-packet techniques. In this chapter, we describe the underlying Hamiltonians and the related reduced dynamical descriptions, show how to recast real-world problems into this form, introduce some of the methods currently used to deal with high-dimensional quantum dynamics, and present the results of this strategy when applied to numerous problems of physicochemical interest.


  • System–bath dynamics
  • Multi-configuration time-dependent methods
  • Generalized Langevin equation
  • Brownian motion
  • Effective modes

1. Introduction

Recent years have witnessed an ever growing interest in dynamical processes that occur in complex environments, for example, ground- and excited-state molecular reaction dynamics in condensed phase, charge and excitation energy transfer in organic functional materials and biomaterials, and elementary processes at the gas–solid interface [16]. Their importance is hardly overemphasized, because of the key role they play in fields as diverse as catalysis, optoelectronics, nanotechnology, biochemistry, and astrophysics, just to mention a few.

The common structure of these problems—a relatively simple system that can be measured and manipulated and that is coupled to an environment only partially under control—has long been subjected to thorough theoretical investigation. Since Einstein’s [7] seminal work on Brownian motion, many important analytical results have been obtained regarding the statistical description of the effect that a large medium—being it a surface or a solvent—has on the dynamics of the small system of interest [8]. The environment, usually designated as the “bath,” is seen to exert two different kinds of force, a friction, and a stochastic force. They result, respectively, in dissipation and fluctuations in the system dynamics and represent two opposite but intimately related effects that ultimately lead to the establishment of equilibrium. This has been made apparent since Langevin formulated the first sound description of an open-system dynamics in 1908, with his equation of motion (the Langevin Equation, LE) and its Generalized version (GLE) [810].

These ideas have been extensively deepened in classical mechanics in the following years and thoroughly validated in both model and realistic systems by several numerical experiments, that is, explicit molecular dynamics (MD) simulations of a “universe” (the system plus the bath) comprising a huge number of degrees of freedom [11].

In quantum dynamics, the situation is considerably more complicated. A brute-force approach is out of reach because of a well-known exponential scaling problem; hence, much of the research in this field focused on reduced dynamical descriptions and aimed at obtaining reasonable master equations for the system density operator [12]. In these open-system quantum dynamical approaches, the degrees of freedom of the environment are traced out and the system undergoes a dissipative, non-unitary dynamics. These approaches are exact in limiting cases only, since most often strong assumptions are needed to obtain manageable equations in a closed form. Among these, the Markov approximation is often invoked (and an effective coarse graining of the dynamics performed) on the basis that the environment correlations last much less than the characteristic time of the system dynamics, a condition that not always holds in practice. Lifting these constraints is possible, for instance, with the help of auxiliary density matrices, but a price of an enormous increase of complexity.

An intermediate possibility between the impractical brute-force approach and the limited reduced dynamics is the so-called system–bath approach, whereby one introduces a surrogate for the environment and explicitly describes its degrees of freedom in the dynamics. In this effective description, the bath is a collection of simple systems (e.g., harmonic oscillators [1315]) and evolves with a relatively simple dynamics. Hence, it is possible to exploit the progress that has been made in the last 20 years or so in propagating high-dimensional wave packets in time [16]. These wave-packet approaches make the unitary evolution of the universe (the system plus the bath) computationally feasible for a large number of coordinates, in many cases large enough to mimic true dissipative environments. The expectation values of interest can then be extracted from the full dynamical evolution using the relevant system operators, and thermal effects can be handled by sampling the mixed initial state of the universe.

In this chapter, we present the work done in the last few years in mapping a physical problem of interest into a system–bath model—the so-called independent oscillator (IO) model—and in solving such model with multi-configuration wave-packet approaches. The dissipative processes investigated range from small amplitude, damped vibrations in model anharmonic systems to “real-world” problems such as hydrogen atom sticking on graphene.

The chapter is organized as follows: Firstly, we describe the IO model in the framework of classical and quantum statistical mechanics, with a focus on its relationship to the generalized Langevin equation and on the role played by the so-called spectral density (SD) of the environmental coupling J0(ω) [13,14,1719]. In addition to the well-established results, we include some recent developments that improved the numerical appeal of the model [2028]. Secondly, we address the problem of mapping a complex (realistic) dynamics into an IO model and deriving the appropriate SD [2931]. We focus in particular on dynamical approaches that suit well to the current practice of accessing dynamical information using on-the-fly simulations, thereby bypassing the need of computing and fitting accurate model potentials. The third part of this contribution deals with the problem of the dynamical evolution of a small system coupled to its IO bath. We describe the now standard and numerically exact multi-configurational time-dependent Hartree (MCTDH) approach [16,32,33], as well as related approximate approaches which better suit to the description of a bath in the IO form. In particular, we discuss in some detail the local coherent-state approximation (LCSA), where the bath evolution is described by a number of Hartree products of pseudo-classically evolving coherent states [3436]. All these approaches are variational and, as such, they share a number of highly desirable features which will be described in some detail. Finally, we present the results of some numerical investigations on both model and realistic systems. Issues, such as vibrational relaxation, decoherence, and scattering, have been extensively investigated in model systems (a harmonic or a Morse oscillator coupled to an oscillator bath) [27,29,3439] and will be summarized in the following. Work on “real-world” problems is still at its infancy, but already offers some notable examples. For instance, we have recently settled a long-standing issue concerning chemisorption of hydrogen atoms on graphene and obtained the first fully quantum and numerically converged results for the probability that the atoms stick on the surface [40,41]. We describe these first exciting results and further provide an outlook of the application of our strategy to other challenging physicochemical problems.


2. Independent oscillator models

The IO Hamiltonian is a popular and extremely powerful tool to study the dynamics of an open system in a quantum setting [1315]. Here, we discuss its connection with the generalized Langevin equation, emphasizing the role played by the SD in the mapping between the two [18,19]. Later, we introduce an effective mode transformation that casts the IO bath into a linear-chain form which suits well to truncation schemes [2022,25,26,28].

2.1. Generalized Langevin equation and SD

The generalized Langevin equation describes the dynamics of a Brownian particle in both the classical and the quantum (Heisenberg) setting. It is a stochastic equation for a system degree of freedom s of mass m subjected to a deterministic potentialV, a random Gaussian force ξ and a friction term, determined by a memory kernel γ(t),

m s ´ t + m γ t τ s . τ d τ + V s t = ξ t E1

Causality of the memory kernel (γ(t) = 0 for t < 0) has important implications for the analytic properties of its Fourier transform γ ~ ω = γ t e i ω t d t when continued to the upper-half complex ω-plane. The SD of the environmental coupling, J0(ω), is defined with the help of the real part of γ ˜ ω

J 0 ω = m ω R γ ω E2

and fully determines the memory kernel by virtue of the Kramers–Kronig relations, namely through

γ t = 2 π m Θ t 0 J 0 ω ω cos ω t d ω E3

where Θ(t) is the Heaviside step function. It further determines the stochastic process ξ(t), provided is Gaussian, by virtue of a fluctuation–dissipation (FD) theorem of the second kind,

Here and in the following  … ⟩denotes an average over the canonical equilibrium.

ξ t ξ 0 = π + J 0 ω 1 e β ω e i ω t d ω E4

here written for a quantum environment (the classical result follows in the limit of high temperatures, β = (kBT)− 1 → 0 and takes the form ⟨ξ(t)ξ(0)⟩ = Θ(t)(t)/β). Hence, all the environment-related terms included in the GLE are uniquely defined by the SD.

Once J0(ω) is known, it can be used to construct an IO Hamiltonian

This is also known as Caldeira–Leggett Hamiltonian, after the seminal work by Caldeira and Leggett on the effects of dissipation on quantum tunneling [17].

H IO = p s 2 2 m + V s + k = 1 p k 2 2 m k + m k ω k 2 2 x k c k m k ω k 2 f s 2 E5

that can be made (quasi) equivalent to the GLE above by appropriately choosing the harmonic oscillator (HO) frequencies and the coupling coefficients. In Eq. (5), the system degree of freedom s is coupled to a collection of harmonic oscillators (xk, pk) of mass

In the following we will adopt, without loosing generality, the same mass for all the oscillators, i.e. mk≡μ for all k, where μ is a numerically convenient value.

mk and frequency ωk. The system–bath coupling is a linear function of the bath coordinates, whereas its dependence on the system coordinate is here specified through the function f(s), which typically complicates the GLE (by introducing state-dependent friction) but is often necessary on physical grounds. The simple linear coupling which makes Eq. (5) equivalent to the standard GLE of Eq. (1) is appropriate to cases where only the near equilibrium configurations of the system are explored (s ≈ 0), but is clearly limited because it describes a coupling which steadily increases when moving the system out of its equilibrium position. Thus, for instance, in previous scattering calculations using a model Morse potential, a coupling function with a finite limit for s → + ∞ [27,36],
f s = 1 e α s α E6

(but yet such that f(s) ≈ s for s ≈ 0 was used to correctly describe the asymptotically free system.

Notice that in Eq. (5), a system potential counter-term f(s)2 appears which balances the distortion induced by the system–bath coupling and ensures the thermodynamic stability of the Hamiltonian [14]. Indeed, in the form given in Eq. (5), the bath adds only quadratic terms (the sum on the r.h.s.) to the system Hamiltonian, thereby guaranteeing a lower bound to the energy spectrum for any reasonable system potential.

The equivalence between the two dynamical formulations [Eqs. (1) vs. (5)] is established when the coupling coefficients sample the SD J0(ω) of the problem, for example, for evenly spaced bath frequencies ωk = kΔω, when the coefficients are set according to

C k = 2 π μ ω k Δ ω J 0 ω k E7

It rigorously holds for a finite time only, determined by the size of the bath in Eq. (5). For longer times, Eq. (1) keeps on describing a dissipative dynamics, whereas the Hamiltonian dynamics of Eq. (5) displays the consequences of discretely sampling the bath frequencies. More precisely, the equivalence is guaranteed up to the Poincaré recurrence time tP = 2π ∕ Δω of the finite system, which thus needs to be set larger than the any time scale of interest. This has to be done by choosing the appropriate discretization Δω, compatibly with a reasonable number of oscillators N and a high frequency cutoff ωc = NΔω of the bath. The latter determines the smaller time that can be resolvedtc = 2π ∕ ωc; higher frequencies, if present, can always be absorbed in a mass renormalization term provided we are not interested in times smaller than tc.

Unlike the starting GLE, the IO Hamiltonian of Eq. (5) can be quantized by applying standard quantization rules. Furthermore, the relatively trivial dynamics followed by the harmonic oscillators of the bath makes the use of standard time-dependent wave-packet approaches to the dynamics possible. As long as a system–bath Hamiltonian can be effectively mapped into a GLE, this represents a powerful and general methodology to tackle an open-system quantum dynamical problem.

2.2. Chain representation of the bath

Recent work has shown that the IO Hamiltonian of Eq. (5) can be expressed in an alternative representation which is particularly suited to truncation schemes or hierarchical descriptions of the dynamics [20,23,24,27,42]. The idea is incredibly simple and reads as follows. Eq. (5) naturally introduces a collective or “effective” mode X 1 = k = 1 N c k x k / D 0 (here D 0 2 = k = 1 N c k 2 is just a kind of normalization constant) that allows one to write the interaction term as

H i n t = K = 1 N c k x k f s = D 0 X 1 f s E8

The definition of this effective mode fixes the first “column” of an orthogonal transformation of the original bath coordinates into a new set of coordinates, otherwise arbitrary. The rest of the transformation matrix can be fixed by requiring that the “residual bath” is in normal form. In this way, Eq. (5) becomes an equivalent IO Hamiltonian for the s plus X1 degree of freedom, coupled to a bath of N − 1 oscillators. The coupling only occurs through X1, and allows one to define a new function J1(ω), that is, the SD felt by the mode X1 as a consequence of its interaction with the residual bath. In the continuum limit, this procedure can be indefinitely iterated to define a sequence of effective modes X1X2 … XM … coupled in a linear chain form and a corresponding sequence of SDs J1(ω), J2(ω) … JM(ω) … which describe the residual bath “felt” by each mode. The sequence of SDs is determined by a simple recurrence relation which can be started with J0(ω) [23,24], that is, Jn + 1(ω) only requires the previous SD Jn(ω) and two of its functionals

+ ω W n J n + 1 ω = D n 2 J n ω E9

Here, + ω W n is the limit on the real axis (from above) of the “Cauchy transform” of Jn(ω),

W n z = 1 π + d ω J n ω ω z E10


D n 2 = 2 π 0 ω J n ω d ω E11

determines the coupling coefficients between the nth and the (n+1)th effective modes. The linear chain form of the Hamiltonian further requires the frequency Ωn + 1 of the (n+1)th effective mode,

Ω n + 1 2 = 2 π D n 2 0 ω 3 J n ω d ω E12

and reads, for mass-scaled coordinates, as

H chain = p s 2 2 m + V s + Δ V s D 0 f s X 1 + n = 1 1 2 P n 2 + 1 2 Ω n 2 X n 2 D n X n X n + 1 E13

Here, ΔV(s) is the potential counter-term

Δ V s = 1 2 m δ Ω 0 2 f s 2 m δ Ω 0 2 = 2 π 0 J 0 ω ω d ω E14

And {(XnPn)}n = 1 … ∞ is the set of effective modes and their conjugated momenta. Equivalently, Eq. (13) can be rearranged to explicitly display its thermodynamic stability [29]

H chain = p s 2 2 m + V s + n = 1 P ' n 2 2 μ + μ ω n 2 2 X ' n β n X ' n 1 2 E15

by introducing new coordinates, X ′ 0 = s and X n = X n / μ for n=1,2.., the bare effective mode frequencies ωn2 = Ωn2 − δΩn2 (where δΩn2 is defined analogously to δΩ02 above, Eq. (14), but in terms of Jn and μ) and the coefficients β 1 = D 0 / μ ω 1 2 and βn = Dn − 1/ωn2 for n > 1.

One interesting issue concerning Eq. (13) is whether a limiting residual SD exists and, in that case, which forms takes lim n J n . It can be readily shown [23] that, if the limit exists, it is the quasi-Ohmic SD given by the Rubin dissipative model. Numerical tests have further shown that this convergence is fast enough that inclusion of a relatively small number of effective modes in an enlarged system makes the resulting dynamics effectively Markovian. Thus, the chain transformation provides a powerful tool for describing non-Markovian phenomena by means of Markovian master equations (applied to enlarged systems).

Strictly speaking such “Markovian reduction” rigorously holds in classical mechanics only; in a quantum setting the very definition of Markovian dynamics is still debated. Thus, one should better refer here to an “Ohmic embedding”.

Further analytical work rigorously proved the convergence conditions of the sequence of spectral densities and established the connection between the chain construction and the moment problem, the theory of orthogonal polynomials, and the Padé approximants [21].

One further issue on the effective mode construction is of much practical interest and concerns its role in defining approximate representations of the bath [24]. In fact, the construction of the linear chain amounts to “unroll” the memory kernel in time—any excitation initially localized in the system necessarily moves sequentially along the chain starting from its end attached to the system. This is contrast to what happens with Eq. (5), where the coupling pattern is appropriate for a frequency resolution of the kernel. As a consequence, truncated or Markov-closed chains with n effective modes can be shown to exactly reproduce γ(t) to the fourth order in time, up to an irrelevant constant of order γ(0)/n [24]. Apart from its conceptual significance, this property proves to be extremely useful for numerical simulations too, since it allows one to single out those bath modes which are most important for the system dynamics. In this representation, in fact, one can easily identify “primary modes” that need to be treated at a high correlation level and “secondary modes” that can be left weakly correlated or even uncorrelated. As it will be shown in the following with some numerical examples, this simple prescription offers the opportunity of tackling long-time issues in explicit dynamical studies of system–bath problems [27].


3. Mapping of a complex system to an IO model

The possibility of using the IO Hamiltonian for simulations of “real-world” systems relies on the existence of a general strategy to derive a GLE from a given microscopic model. In the past, the fundamental problem of mapping in an exact way a reduced dynamics into a GLE was addressed by many authors,

The problem is essentially classical in nature, since the statistical properties of the bath (when subsumed in the spectral density J0(ω) are the same for both the classical and quantum GLE.

but seldom checked in realistic physicochemical problems [43,44]. A substantial contribution in this direction has been recently given by Ivanov and co-workers [30,31], who proposed a classical MD-based methodology that makes use of the combined information of two correlation functions to extract the SD of the bath, and validated it in realistic molecular problems. Their technique is similar in spirit to the one that will be described in the following and that we independently devised at the same time. Though of more limited applicability, the latter has the advantage on relying only on simple dynamical information and does not require any a priori knowledge of the system potential. It thus suits well to an ab initio MD determination of the environmental forces. This inversion procedure (from the dynamics to the SD that generated it) is briefly described in this section, along with an account of its numerical performance [29]. In concluding this section, we further address the problem of going beyond the simple IO model in order to capture anharmonic effects of the environment [29,39].

3.1. Inversion of classical dynamics

When the Brownian motion along s is limited to near-equilibrium configurations, the bare system potential is harmonic, V s = 1 2 m ω 0 2 s 2 , and the SD determines not only the correlation function of the environmental fluctuations [see Eq. (4)], but also the frequency-dependent autocorrelation function of the position,

C ˜ ω = + e i ω t C t d t = + e i ω t s t s 0 d t E16

namely through

1 2 ω C ˜ ω = k B T m I m 1 ω 0 2 ω 2 i ω γ ˜ ω E17

that follows from Eq. (1) by harmonic analysis. This equation relates the reduced dynamics of the system to its coupling with the environment and can be “inverted” to give an analytic expression for J0(ω) in terms of C ˜ ω . This can be accomplished by introducing the retarded correlation function + t = Θ t C t C and by exploiting the analytic properties of its Fourier transform (see Ref. [29] for the detailed derivation). The resulting formula for J0(ω) is as follows

J 0 ω = k B T 2 ω C ˜ ω | Γ + ω | 2 E18

Where ϵ 0 + Γ ω + i ϵ + ω = lim Γ is the limit on the real axis from the u.h.p. of the function,

Γ z = 1 π + ω C ˜ ω / 2 ω z d ω E19

that is, the Cauchy transform of the function f ω = ω C ˜ ω / 2 . Here, the displacement autocorrelation function C(t)(or, equivalently, the velocity autocorrelation) ⟨(t)(0)⟩ = − d2C(t)/dt2 is readily available from equilibrium MD and is the only dynamical information required.

3.2. Numerical tests

In Ref. [29], we tested the inversion procedure of Eq. (18) on some model systems to elucidate how its effectiveness is affected by the presence of an harmonicities and/or by a “Debye” cutoff frequency in the environment. A number of IO Hamiltonians of the kind of Eq. (5) and its variants were used to generate the dynamics, with a reasonable choice of the parameters that mimicked typical molecular problems. Some key alternatives were considered (e.g., a harmonic vs. an anharmonic system oscillator, Ohmic vs. non-Ohmic baths, small vs. large oscillator frequency [compared to the Debye cutoff], etc.), and the position correlation function C(t) was computed with MD and then inverted according to Eq. (18). Some results of these numerical tests are displayed in Figure 1.

The main conclusions of such numerical analysis can be summarized as follows. When the system is a harmonic oscillator, the transformation perfectly recovers the original SD up to the bath Debye frequency ωD (1000 cm− 1 in our numerical tests). For higher frequencies, the spectrum is not identically zero but rather shows an increasing baseline due to the numerical implementation of the Cauchy transform

Numerical evaluation of Eq. (19) requires the introduction of a high-frequency cutoff. The problem arises when using an “unbiased” cut-off frequency well above the spectral range of interest (ωc = 4000cm− 1 in the simulations) and can be easily amended by setting ωc equal to the bath Debye frequency (if known).

of Eq. (19). Further, when the system frequency ωs lies above ωD, the SD features a single Lorentzian peak, which stands out from the background and is robust against variations of the bath and/or the temperature. This is the numerical representation of a Dirac-δ contribution

A fictitious coupling to the bath appears here because numerically the autocorrelation function needs to be damped.

that follows from Eq. (18) when ωs lies outside the support of J0 (ω). Anharmonicity has two strong effects in this case: on the one hand, it causes the appearance of higher harmonics at a frequency 2ωs; on the other hand, it induces a broadening of the Dirac-δ signal that is about twice as large as the spectral width of the bath (i.e., the appearance of combination frequencies ωs ± ωB). Importantly, the temperature has a dramatic effect on the computed SD: At “high” temperature, and in contrast to the HO case, the remnant of the δ-peak undergoes substantial broadening, in striking contrast with the expectation of a temperature-independent SD.

Figure 1.

The spectral densities (atomic units) obtained by “inverting” the dynamical information of the position correlation function according to Eq. (18) (color coded as in that figure) are compared to the original non-Markovian SD used to define the models (green lines). Results are shown for different temperature—50 and 300 K—and for different choices of the system potential—harmonic and Morse—with an intrinsic frequency either above or below the Debye frequency ωD of the bath. Dotted lines mark the frequencies of the system oscillators.

The presence of these features warns against blind application of the transformation of Eq. (18), particularly when the system frequency is larger than ωD. The temperature-dependent background in the SD at high frequencies is unphysical and reflects just the anharmonicity in the system. However, provided such effect is clearly identified, no real problem arises in modeling since the anharmonicity in the system potential can always be included in the IO Hamiltonian by using the appropriate potential V(s).

3.3. Nonlinear extensions of the IO Hamiltonian

As mentioned in the previous section, system anharmonicity poses no real problem to modeling, and only generates a spurious high-frequency coupling that can be minimized by working at low enough temperature. In realistic situations, however, structured features in the spectral region ω > ωD are expected quite generally also from the failure of the bilinear coupling model and/or of the harmonic approximation of the bath oscillators. In both cases, the apparent coupling at frequencies above the Debye bath cutoff does have a physical origin, and thus the question arises whether J0(ω) at such frequencies can be a surrogate for a more complicated coupling model. This is a rather intricate issue related to the general problem of whether a mapping of the system dynamics to a GLE exists and how it can be realized in practice (see, e.g., Ref. [30] for a recent, in-depth analysis of this issue). In general, real molecular oscillators are not harmonic and the system–bath coupling is not bilinear—especially if highly excited vibrational states are being probed—two factors which are hard to disentangle at finite temperature. In these situations, the present “dynamical” approach, when considered in the low-T limit above, can only provide the small-amplitude expansion of the coupling term and needs to be integrated with some empirical knowledge about the interaction between the system and the environment that could guide the formulation of an extended IO model.

In Eq. (5), our definition of IO Hamiltonian, we have already included a coupling which is not linear with respect to the system coordinate. It involves a shape function f(s) that can be used to modulate the strength of the coupling to the bath in the state space of the system, thereby giving rise, as already mentioned, to state-dependent friction. Such extension seems to be a necessary (and simple) modification to address realistic situations; for instance, if s is the height of an adsorbate above a surface, the coupling should vanish for large s and exponentially increase for smalls, as indeed occurs with Eq. (6).

As for the bath, the consideration of a nonlinear coupling poses more problems. In general, an exponential interaction model seems to be appropriate in typical physicochemical situations, where relaxation occurs as a consequence of close encounters between the molecular system of interest and the atoms of the environment. One simple ansatz of this kind is the replacement of f(s)∑ckxk in Eq. (5) with

μ D 0 e α X ' 1 e α X ' 1 1 α f s E20

Such coupling can be justified in the context of the linear-chain representation of the bath (see Section 2.2) and follows from Eq. (15) upon replacing the first harmonic term of the series with a Morse potential with the same frequency. Here, α− 1 is an empirical parameter that describes the characteristic length of the interaction and, for consistency, should be of the order of the atom dimensions. The thus-defined exponential model makes use of the spectral properties of the proper bath only (i.e., for frequencies ω < ωD) to introduce multiphonon relaxing pathways already at the lowest order in perturbation theory, as a simple calculation of the Fermi’s golden-rule rate shows. It is further simple enough to be easily handled with the same numerical methods to be described below and can thus be considered beyond the limits of applicability of perturbation theory.

As a last possibility, there may be realistic situations in which the dynamics of the bath per se shows strong anharmonic effects. For such purpose, in Ref. [39], a new IO Hamiltonian was proposed in which a bath of Morse oscillators was nonlinearly coupled to the system, according to

H Morse = p s 2 2 + V s + k = 1 N p k 2 2 + D k 1 e α k x k c k 2 D k α k s 2 D k E21

where, as before, (sps) and (xkpk) denote the coordinates and momenta of the system and of the bath degrees of freedom, respectively. In Eq. (20), the parameters of the Morse potentials, Dk and αk, can be specified by fixing the oscillator frequencies and ensuring a uniform number of bound states Λ over the bath oscillators. Similarly, the coupling coefficients can be chosen to sample a given SD according to Eq. (7).

The bath dynamics entailed by Eq. (21) can be highly nonlinear, especially for those low-frequency oscillators which undergo large amplitude motion.

In principle, such model also describes energetic processes that irreversibly modify the environment, a phenomenon that can be mimicked by the dissociation of one or more oscillators.

In Ref. [39], this Hamiltonian was used to show that the anharmonicity of the bath induces in fact nontrivial variations in the (quantum) vibrational dynamics and in the corresponding relaxation rates.


4. Techniques for high-dimensional wave-packet dynamics

As is well known, the brute-force numerical solution of the Schrodinger equation rapidly runs into troubles when increasing the system dimensions. This is clearly seen by a simple scaling argument: With N degrees of freedom, if each of them requires on average a basis functions (or grid points), the total number of vector components is aN and that of any operator (matrix) is aN × aN. For a = 6 and double precision (complex) arithmetic, this means ∼ 102N byte for just storing the vector representing the wave function, an incredibly large number for all but the smallest N. Hence, it became clear soon that developing approximate methods was a necessity for facing large-dimensional quantum problems, and much effort has been spent to this end since the dawn of molecular quantum dynamics. The ideal method should be accurate enough to provide reliable information on the dynamics but also sufficiently cheap to be used for extensive sampling of “initial conditions.”

Clearly, one central point to address when devising a method that could fulfill the above requirements is how to define approximate equations of motion. According to our experience, the most general strategy relies on the use of a variational principle that, apart from being physically transparent and mathematically sound, endows the resulting scheme with nontrivial properties. These features have important consequences in practice (e.g., norm and energy conservation issues are settled at the outset) and will be discussed in some detail in the next subsection. In the rest of the section, a number of multi-configurational methods of increasing simplicity (and decreasing computational costs) will be introduced, with a focus on the dynamics at T = 0K where a wave-function suffices. Finite temperature situations can be handled (at least in principle) by applying the same methodologies to the realizations used to sample the mixed (initial) state of the whole system.

4.1. Variational principle and Hamiltonian flows

  1. The time-dependent variational principle is usually stated in terms of the Dirac–Frenkel condition

δ Ψ | i t H | Ψ = 0 E22

where Ψ⟩ is the approximate solution sought for at a time t (the “trial” wave function which lies in some specified variational manifold) and δΨ is an arbitrary variation (an arbitrary vector tangent to the manifold in Ψ⟩). The condition does lead to the Schrödinger equation when the manifold extends over the whole Hilbert state of the system

In fact, this is the condicio sine qua non for the existence of a variational principle.

and admits a very simple interpretation, which is the best seen when multiplying Eq. (22) by –i/. In fact, − i/ℏHΨ⟩ is easily recognized to be the exact time derivative for the system state Ψ⟩, and thus Eq. (22) is seen to be a condition on the “error” vector t - t exact Ψ . As such, it has a trivial geometrical solution: ∂tΨ⟩ needs to be the projection of the exact time derivative onto the tangent space, a recipe that guarantees that the state vector keeps on staying on the manifold during the time evolution.

The above analysis shows that the time-dependent variational principle provides a local-in-time approximation to the system dynamics: It represents the best one can do in the short run, for the given state at time t. It does not offer, though, any guarantee that the solution Ψ(t)⟩ at a finite time t is “close enough” to the exact solution for some specified initial state Ψ(0)⟩, not even that is the best one can do with the specified manifold at that time.

In fact, the best approximation would just be the point on the manifold that lies closest to the exact solution at time t (whose identity may further depend on the adopted metrics).

However, the expectation is that if the trial manifold is flexible enough to include much of the state space spanned by the exact solution during its journey, the variational solution will remain a rather good approximation at any time, and one can thus exploit the nice features associated with the variational principle.

Energy and norm conservation follow immediately from Eq. (22) under quite mild conditions. For the energy, just notice that ∂tΨ is always an allowed variation (and assume H be self-adjoint, of course). Norm conservation follows similarly when the manifold is scale-invariant, since in that case, Ψ itself is an allowed variation. A much deeper result emerges though when the Dirac–Frenkel condition is replaced by the (quasi) equivalent lowest action principle

This is the true variational principle, i.e. a stationarity condition of some cleverly designed functional (something which Eq. (22) is not). However, the two formulations can be shown to be equivalent under quite mild conditions that are usually satisfied in practice.

δS = 0 and can be formulated as follows: Any variational quantum method, under quite mild regularity conditions, can be recast in the form of a symplectomorphism on a symplectic manifold [45]. Here, the relevant action functional is determined by the (real) Lagrangian
L = i 2 Ψ | t Ψ t Ψ ´ | Ψ Ψ | Ψ Ψ | H | Ψ Ψ | Ψ E23

We expand somehow on this issue in the following since it does not seem to be widely appreciated, despite its deep significance and its potential practical utility. To this end, let us first briefly introduce the concept of Hamiltonian flows and symplectic manifolds [46]. A symplectic manifold is a differentiable manifold equipped with a closed, non-degenerate 2-form ω. In a coordinate system xi, it can be written asω = ∑i,j >iωijdxidxj, where dxi are the fundamental 1-form—that is dxi(v) = vi for any tangent vector v in a given point x—and the product of differentials is the so-called “wedge” product. Non-degeneracy means in practice that ωij is non-singular everywhere in the manifold,

This condition restricts the analysis to even dimensional manifolds.

and this allows one to set up a one-to-one map between tangent and co-tangent vectors (1-form). That is, for a given 1-form α = ∑αidxi, there exists an associated vector field Xα such that X α , v α v = ω , and its flow, defined by the curves x ´ i = X α i . Then, given a smooth function H (which can be called a Hamiltonian) and its 1-form dH, the flow induced by its associated vector field XH (what can be called a Hamiltonian flow) conserves the function itself, dH(XH) = ω(XHXH) = 0. Closedness ( = 0 where d is the “exterior” derivative) means that these properties can be “transported” over the manifold and guarantees that the symplectic form ω itself is invariant under any Hamiltonian flow (formally L X H ω = 0 , LY being the Lie derivative along the vector field Y. This forms the basis for Liouville’s theorem).

With this premise in mind, we (smoothly) introduce a set of real variational parameters x, forming a coordinate system in the sample space.

That is, we set Ψ⟩ = Ψ(x)⟩, where x are, e.g., the real and the imaginary parts of a set of complex parameters specifying the wave function. Though not necessarily finite in number or numerable, it is conceptually easier to think of a large but finite number of parameters.

In terms of this parametrization, the Lagrangian reads as
L = i = 1 n x ´ i Z i x H x E24

Where H(x) = ⟨Ψ(x) ∨ H ∨ Ψ(x)⟩/⟨Ψ(x) ∨ Ψ(x)⟩ and Z i = i ħ 2 Ψ Ψ Ψ Ψ x i - Ψ x i Ψ . The latter are the components of a 1-form α = − ∑Zidxi that can be differentiated to give a closed 2-form, ω =  = ∑ωijdxidxj. Provided ωij = ∂Zi/∂xj − ∂Zj/∂xi is non-singular, such 2-form is non-degenerate and thus represents a symplectic structure. In fact, the variational equations of motion take the form

x ´ i ω i j = H x j E25

or, equivalently, for a generic tangent vector v,

x ´ i ω i j v j = ω x ´ v = H x j d x j v = d H v E26

Hence, if ω is a symplectic form, the “variational flow” is the Hamiltonian flow of the Hamiltonian H(x), namely x ´ = X H . Under such circumstances, the equations of motion can also be written with the help of the Poisson brackets

x ˙ i = H x i E27

which are defined by {fg} = ω(XgXf) for any two smooth functions f and g, or, equivalently, by

f g = f x i g x j ξ i j E28

(where ξij is the matrix inverse of ωij) when a coordinate system is introduced.

The importance of the symplectic structure thus described is hardly overemphasized. Apart from its possible consequences on fundamental issues such as the emergence of classicality, or the coexistence of quantum and classical worlds, it offers in practice the possibility of introducing robust propagation schemes in solving the variational equations of motion. These symplectic propagators would not only conserve energy (and norm) but also the whole symplectic structure, a property that might be of great help when numerically investigating the emergence of irreversibility in Hamiltonian systems like the ones described by Eq. (5).

4.2. MCTDH, G-MCDTH, LCSA, and related methods

In this section, we give a brief account of (wave packet) quantum dynamical methods that have been applied in the past to system–bath problems of the kind discussed in this Chapter. The presentation is necessarily limited and, following authors’ personal experience, focuses on the so-called multi-configurational methods only. In these methods, the wave function is written as a superposition of “Hartree products” of single-particle states (or single-particle functions, spf’s), and both the coefficients of such superposition and the single-particle states (or some of them) are variationally optimized. The accuracy, the numerical complexity and the target problems of the method strongly depend on the choice of spf’s, that is, whether they are fully flexible, constrained to a given functional form or frozen. As a result, several different methods exist which stem from the same multi-configurational ansatz.

Among the various possibilities, the conceptually simplest method is the multi-configuration time-dependent Hartree (MTDH), developed decades ago by Meyer et al. [32]. The ansatz is a straightforward expansion of the wave function in terms of (orthonormal) single-particle states,

Ψ = J A J ϕ j 1 1 . . ϕ j k k . . ϕ j F F E29

where the sum runs over the possible configurations labeled by the multi-index J = (j1,.. jN)and ϕ j k for j = 1,.. Nk is a set of orthonormal states for the kth “particle.” Here the “particle” can be either a single degree of freedom or a collection of modes, according to what is called a “mode combination.” The corresponding variational manifold is obtained by varying the complex amplitude coefficients of the superposition and the single-particle states, which are assumed to be orthogonal to each other but otherwise fully flexible in their respective Hilbert space. The orthogonality condition proves to be a key strength of the method, since it guarantees that the configurations entering Eq. (29) are orthogonal to each other.

The equations of motion follow from the application of the Dirac–Frenkel condition to the above MCDTH ansatz, under the orthonormality constraints on the spf’s which are typically introduced by means of arbitrary Hermitian matrices g(k) (one for each mode) that fix the evolution of the spf’s in the “occupied” space. The important equations of motion are for the evolution of the spf’s in the “unoccupied” space and for the evolution of the amplitude coefficients. Their derivation is straightforward, though lengthy, and the result can be summarized as follows. The amplitude coefficients satisfy a kind of matrix form of the Schrödinger equation in the basis of configurations

i A ´ J = L J | H | L A L - i k = 1 F l = 1 N k g j k l k A J k l E30

only corrected for the “gauge” terms arising from the orthonormality constraints. Here, L = (l1,.. lF) is a multi-index analogous to J above, L⟩ and J⟩ are shorthands for the configurations, jk is the kth index of J, and Jk(l) is the same multi-index J with l replacing jk. The orbital equations, on the other hand, are mean-field like

This expression is generally written as an explicit equation for the orbitals, by adding the above mentioned projection onto the occupied space of the spf velocity.

ρ k - 1 j l H lm k | ϕ m k i 1 - P k ϕ ´ j k = 1 - P k l , m = 1 n k E31

and involve the projection Pk onto the space spanned by the spf’s of the kth mode, the mean-field operators H lm k = Ψ l k | H | Ψ m k and the inverse of the density-matrices

Strictly speaking, ρk is the transpose of the 1-particle density for the kth mode, in the spf’s basis.

ρk, ρ k lm = Ψ l k | Ψ m k , here written in terms of hole wave functions, that is, Ψ m k = a m k Ψ , where a m k annihilates the state m of the kth mode. Details on the method and its numerical implementation can be found in the literature [16,32,33].

Some general considerations are in order. The MCDTH method does not solve the exponential scaling problem of quantum dynamics, but considerably alleviates it since replaces a potentially large number of (static) basis functions with a smaller set of “dynamically optimized” elements. As such, it is extremely flexible and allows a systematic search of the convergence of the solution with respect to the length of the expansion of Eq. (29). In fact, provided large enough computational resources are available (how large depends of course on the problema at hand), the MCTDH solution can be made numerically exact.

A second issue concerns the kind of problems MCTDH may handle. The method is “general purpose” and can tackle arbitrary problems, provided the interaction terms between modes can be reasonably described as (sum of) products of terms involving one mode at a time. This is due to the appearance of the mean-field operators above, whose evaluation requires “tracing” over potentially many degrees of freedom. Apart from this, there exists no limitation in the form of the system Hamiltonian and indeed, MCTDH has been applied with success to a very large number of different problems. The application to system–bath problems to be described below represents just one possible problem where the method applies; further applications can be found in the original research papers and in the extensive review literature [16]. Here, we just mention that a user-friendly, highly efficient, general MCDTH code which takes arbitrary Hamiltonians as input is freely available upon request to the author [47].

A second class of multi-configurational methods is represented by the Gaussian-multi-configuration time-dependent Hartree (G-MTDH) developed a while ago by Burghardt et al. [48]. It is still a general-purpose method that can handle different kinds of quantum dynamical problems, and it is obtained from Eq. (22) by fully or partially replacing the flexible spf’s with Gaussians. As a result, the equations of motion for the Gaussians become classical-like with considerable saving of memory and computer time (in fact, one propagates the few parameters needed to define the Gaussians), at the expense of introducing overlap matrix elements between them.

Though the method has several variants depending on the number of Gaussians introduced, it was originally formulated for system–bath-like problems, where one easily identifies primary modes (to be described at the high, fully flexible level) and secondary, less important modes that can be managed with moving Gaussians. In that case, the equations of motion for the amplitude coefficients and for the fully flexible spf’s of the primary modes are similar to those of the MCDTH above, with minor modifications only, whereas a new set of first-order differential equations appear for the Gaussian parameters [48].

Along this line of thought, LCSA was specifically designed as a local coherent-state approximation [34] to the dynamics of system–bath Hamiltonians of the general form

H = H s y s + H e n v . . q k p k . . ; s E32

where Henv is an “environment” Hamiltonian (comprising the coupling with the system) which is supposed to be local in system coordinate(s) s and approximately harmonic in the bath degrees of freedom (..qkpk..). In this case, the shaping of the wave function relies on the fact that (i) the coupling to the bath is local in system coordinates, and (ii) the bath is approximately harmonic. Upon introducing a set of system discrete variable representation (DVR) states

These are highly localized objects in configuration space which underlie the use of any “grid”. See Refs. [49,50] for a formal introduction.

one expands the wave function as

| Ψ = α C α ξ α Φ α E33

where ξα⟩ is a DVR set for the subsystem coordinates, and Φα⟩ are the resulting local bath states (one for each grid point α used to cover the relevant system configuration space). The latter are then written as products of HO coherent states (CSs), that is,

| Φ α = z α 1 z α 2 . . z α F : = Z α E34

and, as a result, the bath dynamics is described by a set of coupled, pseudo-classical trajectories z α k = z α k t , one for each bath degree of freedom k and system grid point α. The system dynamics, on the other hand, is contained in the time evolution of the amplitude coefficients

Rigorously speaking, the system reduced density matrix ρ also requires the overlap between bath state, ραβ = CαCβZβ ∨ Zα⟩ in the underlying DVR.


Equations of motion can be derived from the Dirac–Frenkel condition, using Cα and z α k as dynamical variables [34]. When using conventional phase factors for the CSs, they take the following form. The “system equation” is a kind of Schrödinger–Langevin equation

i C α ´ = β H α β damp C β + v α eff C α E35

in which the elements of the system DVR Hamiltonian are damped by the overlap between bath states, H α β damp = H α β s y s Z α Z β . The local, effective potential, veff = vlmf + vgauge, contains a “local mean-field” potential

k * , z α k , . . ; s α . . z α v α l m f = Z α H e n v s α Z α = H ord e n v E36

(here, H ord e n v is the environment Hamiltonian operator expressed in terms in a k , a k and normally ordered, that is, with all a k ’s on the left of ak’s) and a gauge potential

k * z ´ α k z α I m v α gauge = - i k = 1 F z α k z ´ α k = k = 1 F E37

which can be explicitly written down with the bath equations below. The bath equations are pseudo-classical

k * , . . ; s α . . z α k z α i C α z ´ α k = β H α β damp z β k - z α k C β + C α H ord e n v a k E38

and contain both a classical (local) force (− i/ℏCα times the second term on the r.h.s) and a quantum one (− i/ℏCα times the first term on the r.h.s) coupling the CSs of the same degree of freedom at different grid points. The latter is essential for a quantum, though approximate, description of the bath dynamics. For a detailed derivation of the equations see Ref. [34], and notice that, in general, [aford(aa)] = ∂ford(aa)/∂a.

This concludes the description of the original LCSA method. Several variants are possible (e.g., replacing the DVR states with energy eigenstates or fully flexible system states) and can be found in the original literature [35,36]. Also, the closely related CC-TDSCF method [49] in which the CSs are replaced by fully flexible functions has been shown to provide essentially the same results as LCSA [36], thereby showing the soundness of the CS approximation for the (local) bath dynamics.

In fact, among the features of LCSA, one key strength of the method is that it reduces the bath dynamics to classical-like evolution, with a number of trajectories that scales linearly with the bath dimensions. This means that the method itself has a power-low scaling with such dimensions, the exponent of this scaling depending (eventually) on the interaction between bath modes. For bath modes coupled to the system only [as in Eq. (5)], linear scaling has been observed and model simulations with tens of thousands of bath degrees of freedom performed on modest computers. This good scaling is in common with mixed quantum-classical methods, which, however, fail to correctly represent the system–bath correlations.

Coupled trajectories also arise in a number of closely related approaches, namely the coupled coherent-state method of Shalashilin and Child [50,51] and the G-MCTDH method mentioned above. The latter, in fact, is strongly connected with LCSA and reduces to it as a limiting case (see Appendix B of Ref. [34]). The main difference between the two is that in LCSA all the configurations are orthogonal to each other, as a consequence of the presence of a different DVR state in each of them. This leads to considerable simplifications in the resulting equations, at the price of a reduced accuracy.

Finally, one interesting property about the pseudo-classical description of the bath degrees of freedom is that it suits well to induce dissipative dynamics into the total system. This can be accomplished by adding a suitably designed friction coefficient η to the bath equations, mimicking the presence of a secondary (infinite, memory-less) bath. More formally, it can be shown that applying the LCSA approximation to a system+bath+secondary bath configuration, a classical approximation to the secondary bath dynamics, and standard assumptions (Ohmic bath in the continuum limit) a friction coefficient appears in the LCSA equations for the system+bath degrees of freedom

The same applies to finite temperature cases where, as expected, both a friction and a fluctuating term appear in the LCSA equations of each realization.

(see Appendix A of Ref. [34]). This possibility has been exploited, especially in conjunction with the need of removing numerical instabilities of the method without altering the system dynamics.


5. Applications

Here, we present some numerical applications of the IO model, starting from simple simulations of a Brownian anharmonic oscillator—used as a testing ground for new dynamical methodologies [34,3638] and/or for different representation of the Hamiltonian [27]. The last part of the contribution will be devoted to a “real-world” application, namely the hydrogen atom dynamics on the graphene surface [40,41].

5.1. Model systems

We consider here a model Hamiltonian describing an anharmonic (Morse) oscillator coupled to a heat bath. A typical problem considered in this context is the small amplitude, damped motion of the oscillator. The initial state is taken in product form, with the bath in its ground state (to mimic relaxation at T = 0K), and the system slightly displaced from its equilibrium position.

Figure 2.

The small amplitude relaxation problem described in the text, for an Ohmic bath with relaxation time γ− 1 = 850, 200 fs (panels from left to right). The energy of the system is computed with standard LCSA (blue line), LCSA coupled to a secondary bath with η− 1 = 12 fs (green), and eLCSA (red), and compared with the MCTDH benchmark (black line).

This type of simulations was used in Ref. [34,36] to test the performances of different quantum dynamical approximations (the LCSA and its energy-local version, eLCSA). In Figure 2, the results of these techniques are shown along with benchmark MCTDH results, for different Ohmic spectral densities sampled with a bath of 80 oscillators. A Markovian exponential decay of the energy was found for all but the strongest coupling case, where some energy oscillations are clearly evident. The graph illustrates the main problems of standard LCSA, an inherent numerical instability related to saturation of the bath. These problems are solved in either its “damped” version [34] or with eLCSA [36]. The good agreement between LCSA and MCTDH is impressive, especially in light of the timing of the calculations (for LCSA only 2–3 min on a standard PC).

Similarly, in Ref. [27], the small amplitude relaxation of the Morse oscillator was used to illustrate the advantages of the chain transformation (Section 2.2). Here, MCTDH was used and different degrees of correlation were introduced along the chain, namely a small number of oscillators were described by a full, many particle expansion, whereas the rest of the chain was treated with one spf per mode. In this way, we exploited the strengths of the linear-chain representation of the bath to enlarge the physical time window of the simulation (i.e., to increase the recurrence time) at a computational cost which scaled only linearly with the chain length [27].

Figure 3.

Small amplitude motion for the two non-Ohmic models of Ref. [27]. Computed system energy is shown for chain baths of increasing level of correlation, specified by Np, that is, the number of fully correlated effective modes of the chain (red, blue, green lines for Np = 5, 10 and 15, respectively). Also shown for comparison the benchmark obtained with the bath in normal form, as in Eq. (5) (black line).

Some results for small amplitude relaxation with the bath in linear-chain form are reported in Figure 3, for different structured SDs. The agreement with the benchmark is rather satisfactory and, as expected, the minor discrepancies were removed by increasing the correlation level. This is true both for the average system energy and for more detailed quantities like the position correlation functions.

The Morse oscillator was also used to model a dissipative scattering event, one in which the system is initially asymptotically free and moves toward the potential well where energy exchange with the bath occurs. Typically, in the interaction region, the wave-packet splits into two parts: One gets trapped in the well and fully relaxes on the long run, while the other returns to the asymptotic region. The first fraction, possibly resolved over the collision energy of the incoming wave packet, defines the “sticking” probability (having in minds problems where the bath represents a surface and a projectile sticks to it).

Figure 4.

Sticking probability as a function of the incident energy, for the Morse scattering problem described in the main text. Left panel: Results for eLCSA (red circles) are compared to the MCTDH benchmark (black circles) for an Ohmic model with γ− 1 = 1 ps. Right panel: The structured SD of Ref. [27] is adopted, and results are shown for baths in chain form with three different correlation schemes of increasing Np, that is, the number of fully correlated chain oscillators (red, blue, green circles for Np = 5, 10 and 15, respectively). Black circles are benchmark results obtained with the bath in normal form. Lines serve as a guide to the eye.

The sticking problem was considered as a test case for both LCSA and for the linear-chain representation of the bath with MCTDH (Figure 3). Simulations with standard LCSA showed that the numerical instabilities were too severe to extract meaningful sticking coefficients, even if the energy dissipation was described quite accurately [35]. On the contrary, the energy-local variant eLCSA gave stable results but only in semi-quantitative agreement with the benchmark. A detailed analysis showed that this is due to an inadequate system–bath correlation in the adopted ansatz, which is crucial for the energy transfer and hence the sticking process. Excellent results, on the other hand, were obtained by applying MCTDH with a partially correlated linear chain of oscillators [27]. Importantly, the results were shown to steadily and quickly converge toward the benchmark when increasing the level of correlation.

5.2. Hydrogen atom dynamics on graphene

In the last decade, the activated dynamics of hydrogen sticking on graphitic/graphenic surfaces has been one of the most studied gas–surface scattering problems. Despite the apparent simplicity of the system, the presence of both dissipative and quantum features makes it a challenging dynamical problem.

Recently, we have devised a rather elaborate system–bath model to describe hydrogen chemisorption of graphene and used it in a fully quantum study of the sticking dynamics with the MCDTH method. The model consists of an accurate description of the hydrogen atom and its bonding carbon atom (a 4D system), which were then coupled to the graphene sheet described by a phonon bath. It rests on the following, reasonable assumptions: (i) The energy exchange that occurs between the system and the lattice for near equilibrium configurations is representative of energy dissipation; (ii) relaxation proceeds through sequential energy transfer from the hydrogen atom to the carbon atom; (iii) a mapping holds, at least approximately, which relates the classical Hamiltonian dynamics of the interesting C and H atoms to a GLE description. On this basis, the following form was adopted for the Hamiltonian

H = p H 2 2 m H + p C 2 2 m C + V s x H z C + k = 1 F p k 2 2 + ω k 2 2 q k - c k ω k 2 z C - z C eq 2 E39

Here, xH is the position of the H atom, zC the height of the binding C atom above the surface, pH and pC the corresponding momenta, and Vs(xHzC) an appropriate 4D system potential. The frequencies ωk and couplings ck of the IO bath were chosen to sample the spectral density J c ( ω ) ) that describes the coupling of the C atom to the rest of the lattice. The latter was obtained by applying twice the inversion procedure (Section 3.1), using as only input the position correlation function CH(t) describing the equilibrium dynamics of the hydrogen atom [40].

The thus-obtained SD JC(ω) is shown in Figure 5 and presents a clear separation between the low-frequency region (0–900 cm− 1) – associated with the “surface” stretching

The “surface” stretching is one of the normal mode of the 4D system potential Vs(xHzC). This eigen-mode lying at ~460 cm− 1 approximately corresponds to block oscillations of the C-H unit above the surface plane.

—and the high-frequency region of the C–H stretching. Details on how it was extracted and thoroughly checked can be found in the original research paper [40].

Figure 5.

The function J c ( ω ) ), that is, the spectral density “felt” by the carbon in the CH model of Eq. (39). The SD has been obtained from the inversion procedure described in the text and used for the high-dimensional quantum dynamics calculations of Refs. [40,41].

Once the coupling of the C atom with its environment was introduced, the Hamiltonian model of Eq. (39) could be used to investigate the hydrogen atom dynamics. We start here by considering the relaxation problem of the C–H bond. In this case, the system was prepared in an eigenstate of the Vs potential and allowed to relax because of the interaction with the bath (initially in its ground state, to mimic a T = 0 K situation).

Figure 6.

Time evolution of the system energy for the lowest lying vibrational states of the CH potential, computed with MCTDH and here labeled by the three appropriate quantum numbers (νs,CHνs,surfνb) where νs,CH is the C–H stretching, νs,surf is the surface stretching, and νb is the doubly degenerate C–H bending. The panel on the left gives the vibrational energy spectrum of the system. The dashed lines mark the recurrence time of the bath models adopted.

As is shown in Figure 5, relaxation from the surface stretching mode proceeds over a very short time scale and is complete in a few tens of fs. Despite the fast relaxation dynamics indicates a strong coupling between this coordinate and the bath, the energy decay shows essentially a Markovian behavior, except for the slippage at short time which extends for a considerable fraction of the relaxation window. This feature is related to the prepared initial states and to the switching on of the coupling term, which actually causes a slight increase of the system energy. The opposite behavior, on the other hand, is seen in the relaxation of the C–H stretching mode that takes place over a picosecond scale and seems to be complete on a time scale much larger than the 3.0 ps limit imposed by the recurrence time of our bath discretization. The resulting relaxation rate (τ ∼ 5.0 ps) is determined only by the background around the main peak in the SD. Its magnitude is (maybe incidentally) very close to the result obtained by Sakong and Kratzer [52], who applied perturbation theory from first principles and found τ = 5.2 ps.

Next, we consider the quantum simulations of the collinear sticking dynamics [41], that is, the process in which a gas-phase hydrogen atom colliding at normal incidence above a carbon atom exchanges energy with the surface and gets trapped in the chemisorption well. We used the MCTDH method once again and, in addition, classical and quasi-classical methods to single out quantum effects in the results (Figure 7).

Figure 7.

Sticking probability as a function of the collision energy as obtained from the MCTDH calculations with the CH system–bath model described in the main text (green curve). Also shown for comparison the classical sticking curves at T = 50 and 300 K (solid and dashed black lines, respectively) and the quasi-classical results (red line), obtained with the same IO model.

Before discussing the quantum dynamical results, it is instructive to first focus on some classical aspects of the sticking process. In Figure 7, the classical results are reported for two different surface temperatures, T =50 and 300 K. At low temperature, the sticking probability is negligibly small below the (static) barrier energy of ∼ 0.24 eV, whereas above the barrier it reaches a saturation value Ps ≈ 1.0 in a relatively narrow energy range, and decreases afterward. The detail analysis of the dynamical process [41] shows that in the below-barrier energy regime sticking is only determined by the probability that the projectile overcomes the barrier, since in the interaction region, the atom easily dissipates the (small) amount of energy required to get trapped in the chemisorption well. Above the barrier, energy transfer to the surface represents the limiting factor to sticking: Only if a large amount of energy can be transferred to the bath, the projectile is prevented to recross the barrier and to return to the gas phase. As a consequence, a simple, impulsive model of the dynamics describes the results rather well [41].

The quantum results differ from the classical ones in the whole energy range considered. While quantum effects (tunneling) can be invoked in the low-energy regime, above saturation the discrepancy is necessarily due to the quantum nature of the low-temperature surface which, in this T = 0 K limit, shows pronounced zero-point energy effects on the projectile dynamics. As a consequence, quasi-classical simulations show a rather good agreement with the quantum results, apart from the threshold region where tunneling through the barrier dominates. A logarithmic plot of the curves in this energy region (Figure 7) shows though that the effect of tunneling is moderate (less than one order of magnitude), in contrast with the effect of the lattice quantum fluctuations (up two orders of magnitudes). Furthermore, it has been shown that the quantum results are well described by the above-mentioned impulsive model of the dynamics, provided it is extended in order to account for the lattice quantum fluctuations and it is applied at energies not too small compared to the barrier height [41].



This work has arisen over the years through an intense and fruitful collaboration with several people. Among them, we are particularly in debt to Mathias Nest, Peter Saalfrank, Keith Hughes, and Irene Burghardt, who are sincerely acknowledged for their contribution. Support for the numerics from the Regione Lombardia and the CINECA High Performance Computing Center through the LISA initiative is also gratefully acknowledged.


  1. 1. Nitzan A. Chemical dynamics in condensed phases—relaxation, transfer and reactions in condensed molecular systems. Oxford University Press; 2006.
  2. 2. Lee H, Cheng Y-C, Fleming GR. Coherence dynamics in photosynthesis: protein protection of excitonic coherence. Science 2007;316:1462–5. doi:10.1126/science.1142188.
  3. 3. Collini E, Scholes GD. Coherent intrachain energy migration in a conjugated polymer at room temperature. Science 2009;323:369–73. doi:10.1126/science.1164016.
  4. 4. Hwang I, Scholes GD. Electronic energy transfer and quantum-coherence in π-conjugated polymers. Chem Mater 2011;23:610–20. doi:10.1021/cm102360x.
  5. 5. Engel GS, Calhoun TR, Read EL, Ahn T-K, Mancal T, Cheng Y-C, et al. Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems. Nature 2007;446:782–6. doi:10.1038/nature05678.
  6. 6. Collini E, Wong CY, Wilk KE, Curmi PMG, Brumer P, Scholes GD. Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature. Nature 2010;463:644–7. doi:10.1038/nature08811.
  7. 7. Einstein A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann Phys 1905;322:549–60. doi:10.1002/andp.19053220806.
  8. 8. Kubo R, Toda M, Hashitsume N. Statistical physics II—nonequilibrium statistical mechanics. Springer Verlag; 1991.
  9. 9. Lemons DS, Gythiel A. Paul Langevin’s 1908 paper “On the Theory of Brownian Motion” [“Sur la théorie du mouvement brownien,” C. R. Acad. Sci. (Paris) 146, 530–533 (1908)]. Am J Phys 1997;65:1079–81. doi:10.1119/1.18725.
  10. 10. Zwanzig R. Nonequilibrium statistical mechanics. USA: Oxford University Press; 2001.
  11. 11. Tuckerman M. Statistical mechanics: theory and molecular simulation. OUP Oxford; 2010.
  12. 12. Breuer HP, Petruccione F. The theory of open quantum systems. Oxford University Press; 2007.
  13. 13.

    Caldeira AO, Leggett AJ. Influence of damping on quantum interference: an exactly soluble model. Phys Rev A 1985;31:1059–66. doi:10.1103/PhysRevA.31.1059.

  14. 14. Ford GW, Lewis JT, O’Connell RF. Quantum Langevin equation. Phys Rev A 1988;37:4419–28. doi:10.1103/PhysRevA.37.4419.
  15. 15. Zwanzig R. Nonlinear generalized Langevin equations. J Stat Phys 1973;9:215–20.
  16. 16. Meyer H-D, Gatti F, Worth GA, editors. Multidimensional quantum dynamics: MCTDH theory and applications. Weinheim: Wiley-VCH Verlag GmbH & Co. KGaA; 2009.
  17. 17. Caldeira AO, Leggett AJ. Influence of dissipation on quantum tunneling in macroscopic systems. Phys Rev Lett 1981;46:211–4. doi:10.1103/PhysRevLett.46.211.
  18. 18. Weiss U. Quantum dissipative systems. 3rd ed. Singapore: World Scientific; 2008.
  19. 19. Pottier N. Nonequilibrium statistical physics: linear irreversible processes. Oxford University Press; 2010.
  20. 20. Hughes KH, Christ CD, Burghardt I. Effective-mode representation of non-Markovian dynamics: a hierarchical approximation of the spectral density. I. Application to single surface dynamics. J Chem Phys 2009;131:024109. doi:10.1063/1.3159671.
  21. 21. Chin AW, Rivas Á, Huelga SF, Plenio MB. Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials. J Math Phys 2010;51:092109. doi:10.1063/1.3490188.
  22. 22. Prior J, Chin AW, Huelga SF, Plenio MB. Efficient simulation of strong system-environment interactions. Phys Rev Lett 2010;105:050404. doi:10.1103/PhysRevLett.105.050404.
  23. 23. Martinazzo R, Vacchini B, Hughes KH, Burghardt I. Communication: Universal Markovian reduction of Brownian particle dynamics. J Chem Phys 2011;134:011101. doi:10.1063/1.3532408.
  24. 24. Martinazzo R, Hughes KH, Burghardt I. Unraveling a Brownian particle’s memory with effective mode chains. Phys Rev E 2011;84:030102(R). doi:10.1103/PhysRevE.84.030102.
  25. 25. Burghardt I, Hughes KH, Martinazzo R, Tamura H, Gindensperger E, Köppel H, et al. Conical intersections coupled to an environment. In: Domcke W, Yarkony DR, Köppel H, editors. Conical Intersect., vol. 17, World Scientific; 2011, pp. 301–46.
  26. 26. Martinazzo R, Hughes KH, Burghardt I. Hierarchical effective-mode approach for extended molecular systems. In: Hoggan EP, Brändas JE, Maruani J, Piecuch P, Delgado-Barrio G, editors. Adv. Theory Quantum Syst. Chem. Phys., Dordrecht, Netherlands: Springer; 2012, p. 269–83.
  27. 27. Bonfanti M, Tantardini GF, Hughes KH, Martinazzo R, Burghardt I. Compact MCTDH wave functions for high-dimensional system–bath quantum dynamics. J Phys Chem A 2012;116:11406–13. doi:10.1021/jp3064504.
  28. 28. Woods MP, Groux R, Chin AW, Huelga SF, Plenio MB. Mappings of open quantum systems onto chain representations and Markovian embeddings. J Math Phys 2014;55:032101. doi:10.1063/1.4866769.
  29. 29. Bonfanti M, Hughes KH, Burghardt I, Martinazzo R. Vibrational relaxation and decoherence in structured environments: a numerical investigation. Ann Phys 2015;527:556–69. doi:10.1002/andp.201500144.
  30. 30. Gottwald F, Ivanov SD, Kühn O. Applicability of the Caldeira–Leggett model to vibrational spectroscopy in solution. J Phys Chem Lett 2015;6:2722–7. doi:10.1021/acs.jpclett.5b00718.
  31. 31. Gottwald F, Ivanov SD, Kühn O. Vibrational spectroscopy via the Caldeira–Leggett model with anharmonic system potentials. ArXiv E-Prints 2016.
  32. 32. Meyer H-D, Manthe U, Cederbaum LS. The multi-configurational time-dependent Hartree approach. Chem Phys Lett 1990;165:73–8. doi:10.1016/0009-2614(90)87014-I.
  33. 33. Beck MH, Jäckle A, Worth GA, Meyer H-D. The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys Rep 2000;324:1–105. doi:10.1016/S0370-1573(99)00047-2.
  34. 34. Martinazzo R, Nest M, Saalfrank P, Tantardini GF. A local coherent-state approximation to system-bath quantum dynamics. J Chem Phys 2006;125:194102. doi:10.1063/1.2362821.
  35. 35. Martinazzo R, Burghardt I, Martelli F, Nest M. Local coherent-state approximation to system-bath quantum dynamics. In: Shalashilin DV, De Miranda MP, editors. vol. Multidimensional quantum mechanics with trajectories, CCP6; 2009.
  36. 36.

    López-López S, Nest M, Martinazzo R. Generalized CC-TDSCF and LCSA: The system-energy representation. J Chem Phys 2011;134:014102. doi:10.1063/1.3518418.

  37. 37. Nest M, Meyer H-D. Dissipative quantum dynamics of anharmonic oscillators with the multiconfiguration time-dependent Hartree method. J Chem Phys 2003;119:24–33. doi:10.1063/1.1576384.
  38. 38. Burghardt I, Nest M, Worth GA. Multiconfigurational system-bath dynamics using Gaussian wave packets: energy relaxation and decoherence induced by a finite-dimensional bath. J Chem Phys 2003;119:5364–78. doi:10.1063/1.1599275.
  39. 39. López-López S, Martinazzo R, Nest M. Benchmark calculations for dissipative dynamics of a system coupled to an anharmonic bath with the multiconfiguration time-dependent Hartree method. J Chem Phys 2011;134. doi:10.1063/1.3556940.
  40. 40. Bonfanti M, Jackson B, Hughes KH, Burghardt I, Martinazzo R. Quantum dynamics of hydrogen atoms on graphene. I. System-bath modeling. J Chem Phys 2015;143. doi:10.1063/1.4931116.
  41. 41. Bonfanti M, Jackson B, Hughes KH, Burghardt I, Martinazzo R. Quantum dynamics of hydrogen atoms on graphene. II. Sticking. J Chem Phys 2015;143. doi:10.1063/1.4931117.
  42. 42. Hughes KH, Christ CD, Burghardt I. Effective-mode representation of non-Markovian dynamics: a hierarchical approximation of the spectral density. II. Application to environment-induced nonadiabatic dynamics. J Chem Phys 2009;131:124108. doi:10.1063/1.3226343.
  43. 43. Mori H. Transport, collective motion, and Brownian motion. Prog Theor Phys 1965;33:423. doi:10.1143/PTP.33.423.
  44. 44. Zwanzig R. Memory effects in irreversible thermodynamics. Phys Rev 1961;124:983–92. doi:10.1103/PhysRev.124.983.
  45. 45. Kramer P, Saraceno M, editors. Geometry of the time-dependent variational principle in quantum mechanics. Lecture notes in physics. vol. 140. Berlin: Springer Verlag; 1981.
  46. 46. Flanders H. Differential forms with applications to the physical sciences. Dover Publications; 1963.
  47. 47. Worth GA, Beck MH, Jäckle A, Meyer H-D. The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007), Version 8.5 (2011). Current version: 8.5.3 (2013). See University of Heidelberg, Germany; n.d.
  48. 48. Burghardt I, Meyer H-D, Cederbaum LS. Approaches to the approximate treatment of complex molecular systems by the multiconfiguration time-dependent Hartree method. J Chem Phys 1999;111:2927–39. doi:
  49. 49. Zhang DH, Bao W, Yang M. Continuous configuration time-dependent self-consistent field method for polyatomic quantum dynamical problems. J Chem Phys 2005;122:091101. doi:10.1063/1.1869496.
  50. 50. Shalashilin DV, Child MS. Time dependent quantum propagation in phase space. J Chem Phys 2000;113:10028.
  51. 51. Shalashilin DV, Child MS. The phase space CCS approach to quantum and semiclassical molecular dynamics for high-dimensional systems. Chem Phys 2004;304:103–20. doi:10.1016/j.chemphys.2004.06.013.
  52. 52. Sakong S, Kratzer P. Hydrogen vibrational modes on graphene and relaxation of the C–H stretch exmixed-citation from first-principles calculations. J Chem Phys 2010;133:054505. doi:10.1063/1.3474806.


  • Here and in the following  … 〉denotes an average over the canonical equilibrium.
  • This is also known as Caldeira–Leggett Hamiltonian, after the seminal work by Caldeira and Leggett on the effects of dissipation on quantum tunneling [17].
  • In the following we will adopt, without loosing generality, the same mass for all the oscillators, i.e. mk≡μ for all k, where μ is a numerically convenient value.
  • Strictly speaking such “Markovian reduction” rigorously holds in classical mechanics only; in a quantum setting the very definition of Markovian dynamics is still debated. Thus, one should better refer here to an “Ohmic embedding”.
  • The problem is essentially classical in nature, since the statistical properties of the bath (when subsumed in the spectral density J0(ω) are the same for both the classical and quantum GLE.
  • Numerical evaluation of Eq. (19) requires the introduction of a high-frequency cutoff. The problem arises when using an “unbiased” cut-off frequency well above the spectral range of interest (ωc = 4000cm− 1 in the simulations) and can be easily amended by setting ωc equal to the bath Debye frequency (if known).
  • A fictitious coupling to the bath appears here because numerically the autocorrelation function needs to be damped.
  • In principle, such model also describes energetic processes that irreversibly modify the environment, a phenomenon that can be mimicked by the dissociation of one or more oscillators.
  • In fact, this is the condicio sine qua non for the existence of a variational principle.
  • In fact, the best approximation would just be the point on the manifold that lies closest to the exact solution at time t (whose identity may further depend on the adopted metrics).
  • This is the true variational principle, i.e. a stationarity condition of some cleverly designed functional (something which Eq. (22) is not). However, the two formulations can be shown to be equivalent under quite mild conditions that are usually satisfied in practice.
  • This condition restricts the analysis to even dimensional manifolds.
  • That is, we set Ψ〉 = Ψ(x)〉, where x are, e.g., the real and the imaginary parts of a set of complex parameters specifying the wave function. Though not necessarily finite in number or numerable, it is conceptually easier to think of a large but finite number of parameters.
  • This expression is generally written as an explicit equation for the orbitals, by adding the above mentioned projection onto the occupied space of the spf velocity.
  • Strictly speaking, ρk is the transpose of the 1-particle density for the kth mode, in the spf’s basis.
  • These are highly localized objects in configuration space which underlie the use of any “grid”. See Refs. [49,50] for a formal introduction.
  • Rigorously speaking, the system reduced density matrix ρ also requires the overlap between bath state, ραβ = CαCβ〈Zβ ∨ Zα〉 in the underlying DVR.
  • The same applies to finite temperature cases where, as expected, both a friction and a fluctuating term appear in the LCSA equations of each realization.
  • The “surface” stretching is one of the normal mode of the 4D system potential Vs(xH, zC). This eigen-mode lying at ~460 cm− 1 approximately corresponds to block oscillations of the C-H unit above the surface plane.

Written By

Matteo Bonfanti and Rocco Martinazzo

Submitted: 30 October 2015 Reviewed: 24 February 2016 Published: 24 August 2016