## Abstract

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.

### Keywords

- 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 [1–6]. 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) [8–10].

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 [13–15]) 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 J_{0}(ω) [13,14,17–19]. In addition to the well-established results, we include some recent developments that improved the numerical appeal of the model [20–28]. Secondly, we address the problem of mapping a complex (realistic) dynamics into an IO model and deriving the appropriate SD [29–31]. 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 [34–36]. 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,34–39] 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 [13–15]. 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 [20–22,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 potential*V*, a random Gaussian force *ξ* and a friction term, determined by a memory kernel *γ*(*t*),

Causality of the memory kernel (*γ*(*t*) = 0 for *t* < 0) has important implications for the analytic properties of its Fourier transform _{0}(ω), is defined with the help of the real part of

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

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,[1] -

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

Once J_{0}(ω) is known, it can be used to construct an IO Hamiltonian[1] -

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 (x_{k}, p_{k}) of mass[1] - m_{k} 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],

(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 J_{0}(ω) of the problem, for example, for evenly spaced bath frequencies ω_{k} = kΔω, when the coefficients are set according to

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 t_{P} = 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 resolvedt_{c} = 2π ∕ ω_{c}; higher frequencies, if present, can always be absorbed in a *mass* renormalization term provided we are not interested in times smaller than t_{c}.

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

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 X*_{1} degree of freedom, coupled to a bath of *N* − 1 oscillators. The coupling *only* occurs through *X*_{1}, and allows one to define a new function J_{1}(ω), that is, the SD felt by the mode *X*_{1} 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 *X*_{1}, *X*_{2} … *X*_{M} … coupled in a linear chain form and a corresponding sequence of SDs J_{1}(ω), J_{2}(ω) … J_{M}(ω) … 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 J_{0}(ω) [23,24], that is, J_{n + 1}(ω) only requires the previous SD J_{n}(ω) and two of its functionals

Here,_{n}(ω),

and

determines the coupling coefficients between the *n*^{th} 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,

and reads, for mass-scaled coordinates, as

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

And {(*X*_{n}, *P*_{n})}_{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]

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

One interesting issue concerning Eq. (13) is whether a limiting residual SD exists and, in that case, which forms takes

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,[1] - 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,

namely through

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 J_{0}(ω) in terms of_{0}(ω) is as follows

Where

that is, the Cauchy transform of the function *C*(*t*)(or, equivalently, the velocity autocorrelation) ⟨(t)(0)⟩ = − d^{2}*C*(*t*)/*d*t^{2} 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[1] - 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[1] - that follows from Eq. (18) when *ω*_{s} lies outside the support of J_{0} (ω). 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.

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 *J*_{0}(*ω*) 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 small*s*, 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*)∑*c*_{k}*x*_{k} in Eq. (5) with

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

where, as before, (*s*, *p*_{s}) and (*x*_{k}, *p*_{k}) 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, *D*_{k} 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.[1] - 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 *a*^{N} and that of any operator (matrix) is *a*^{N} × *a*^{N}. For *a* = 6 and double precision (complex) arithmetic, this means ∼ 10^{2N} *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* = 0*K* 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

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

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[1] - 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}*Ψ*⟩ 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.[1] - 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[1] - *δ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

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 *x*^{i}, it can be written as*ω* = ∑_{i,j >i}*ω*_{ij}*dx*^{i}*dx*^{j}, where *dx*^{i} are the fundamental 1-form—that is *dx*^{i}(*v*) = *v*^{i} 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,[1] - 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 *α* = ∑*α*_{i}*dx*^{i}, there exists an associated vector field *X*_{α} such that *flow*, defined by the curves*H* (which can be called a Hamiltonian) and its 1-form *dH*, the flow induced by its associated vector field *X*_{H} (what can be called a Hamiltonian flow) conserves the function itself, *dH*(*X*_{H}) = *ω*(*X*_{H}, *X*_{H}) = 0. Closedness (*dω* = 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*_{Y} 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.[1] - In terms of this parametrization, the Lagrangian reads as

Where *H*(*x*) = ⟨*Ψ*(*x*) ∨ *H* ∨ *Ψ*(*x*)⟩/⟨*Ψ*(*x*) ∨ *Ψ*(*x*)⟩ and *α* = − ∑*Z*_{i}*dx*^{i} that can be differentiated to give a closed 2-form, *ω* = *dα* = ∑*ω*_{ij}*dx*^{i}*dx*^{j}. Provided *ω*_{ij} = ∂*Z*_{i}/∂*x*^{j} − ∂*Z*_{j}/∂*x*^{i} 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

or, equivalently, for a generic tangent vector *v*,

Hence, if *ω* is a symplectic form, the “variational flow” *is* the Hamiltonian flow of the Hamiltonian *H*(*x*), namely

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

(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,

where the sum runs over the possible configurations labeled by the multi-index *J* = (*j*_{1},.. *j*_{N})and *j* = 1,.. *N*_{k} is a set of orthonormal states for the *k*^{th} “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

only corrected for the “*gauge”* terms arising from the orthonormality constraints. Here, *L* = (*l*_{1},.. *l*_{F}) is a multi-index analogous to *J* above, *L*⟩ and *J*⟩ are shorthands for the configurations, *j*_{k} is the *k*^{th} index of *J*, and *J*^{k}(*l*) is the same multi-index *J* with *l* replacing *j*_{k}. The orbital equations, on the other hand, are mean-field like[1] -

and involve the projection *P*_{k} onto the space spanned by the *spf’s* of the *k*^{th} mode, the mean-field operators *ρ*_{k}, *m* of the *k*^{th} 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

where *H*^{env} 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 (..*q*_{k}, *p*_{k}..). In this case, the shaping of the wave function relies on the fact that (i) the coupling to the bath is *loca*l in system coordinates, and (ii) the bath is approximately *harmonic*. Upon introducing a set of system discrete variable representation (DVR) states[1] - one expands the wave function as

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,

and, as a result, the bath dynamics is described by a set of coupled, *pseudo-classical* trajectories *k and* system grid point *α*. The system dynamics, on the other hand, is contained in the time evolution of the amplitude coefficients[1] - *C*_{α}.

Equations of motion can be derived from the Dirac–Frenkel condition, using *C*_{α} and

in which the elements of the system DVR Hamiltonian are *damped* by the overlap between bath states, *v*^{eff} = *v*^{lmf} + *v*^{gauge}, contains a “local mean-field” potential

(here, *a*_{k}’s) and a *gauge* potential

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

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, [*a*, *f*_{ord}(*a*^{†}, *a*)] = ∂*f*_{ord}(*a*^{†}, *a*)/∂*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[1] - (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,36–38] 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* = 0*K*), and the system slightly displaced from its equilibrium position.

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

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

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

Here, **x**_{H} is the position of the H atom, *z*_{C} the height of the binding C atom above the surface, **p**_{H} and *p*_{C} the corresponding momenta, and *V*_{s}(**x**_{H}, *z*_{C}) an appropriate 4D system potential. The frequencies *ω*_{k} and couplings *c*_{k} of the IO bath were chosen to sample the spectral density *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 *C*_{H}(*t*) describing the equilibrium dynamics of the *hydrogen atom* [40].

The thus-obtained SD *J*_{C}(*ω*) is shown in Figure 5 and presents a clear separation between the low-frequency region (0–900 cm^{− 1}) – associated with the “surface” stretching[1] -—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].

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 *V*_{s} potential and allowed to relax because of the interaction with the bath (initially in its ground state, to mimic a *T* = 0 K situation).

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

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 *P*_{s} ≈ 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].

## Acknowledgments

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.

## Notes

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