Relaxation dynamics in quantum dissipative systems: The microscopic effect of intramolecular vibrational energy redistribution L. Uranga-Piña and J. C. Tremblay Citation: The Journal of Chemical Physics 141, 074703 (2014); doi: 10.1063/1.4892376 View online: http://dx.doi.org/10.1063/1.4892376 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/7?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Accurate ab initio and “hybrid” potential energy surfaces, intramolecular vibrational energies, and classical ir spectrum of the water dimer J. Chem. Phys. 130, 144314 (2009); 10.1063/1.3112403 Stochastic Liouville equations for hydrogen-bonding fluctuations and their signatures in two-dimensional vibrational spectroscopy of water J. Chem. Phys. 123, 114504 (2005); 10.1063/1.2008251 Intramolecular vibrational energy redistribution in the highly excited fluoroform molecule: A quantum mechanical study using the multiconfiguration time-dependent Hartree algorithm J. Chem. Phys. 120, 6992 (2004); 10.1063/1.1668639 Quantum mechanical study of intramolecular vibrational energy redistribution in the third CH stretch overtone state in benzene J. Chem. Phys. 111, 5617 (1999); 10.1063/1.479857 The effect of two- and three-body interactions in Ar n CO 2 (n=1,2) on the asymmetric stretching CO 2 coordinate: An ab initio study J. Chem. Phys. 106, 10215 (1997); 10.1063/1.474105

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

THE JOURNAL OF CHEMICAL PHYSICS 141, 074703 (2014)

Relaxation dynamics in quantum dissipative systems: The microscopic effect of intramolecular vibrational energy redistribution L. Uranga-Piña1,2 and J. C. Tremblay2,a) 1 2

Facultad de Física, Universidad de la Habana, San Lázaro y L, Vedado, 10400 Havana, Cuba Institute for Chemistry and Biochemistry, Freie Universität Berlin, Takustr. 3, D-14195 Berlin, Germany

(Received 24 April 2014; accepted 25 July 2014; published online 15 August 2014) We investigate the effect of inter-mode coupling on the vibrational relaxation dynamics of molecules in weak dissipative environments. The simulations are performed within the reduced density matrix formalism in the Markovian regime, assuming a Lindblad form for the system-bath interaction. The prototypical two-dimensional model system representing two CO molecules approaching a Cu(100) surface is adapted from an ab initio potential, while the diatom-diatom vibrational coupling strength is systematically varied. In the weak system-bath coupling limit and at low temperatures, only first order non-adiabatic uni-modal coupling terms contribute to surface-mediated vibrational relaxation. Since dissipative dynamics is non-unitary, the choice of representation will affect the evolution of the reduced density matrix. Two alternative representations for computing the relaxation rates and the associated operators are thus compared: the fully coupled spectral basis, and a factorizable ansatz. The former is well-established and serves as a benchmark for the solution of Liouville-von Neumann equation. In the latter, a contracted grid basis of potential-optimized discrete variable representation is tailored to incorporate most of the inter-mode coupling, while the Lindblad operators are represented as tensor products of one-dimensional operators, for consistency. This procedure results in a marked reduction of the grid size and in a much more advantageous scaling of the computational cost with respect to the increase of the dimensionality of the system. The factorizable method is found to provide an accurate description of the dissipative quantum dynamics of the model system, specifically of the time evolution of the state populations and of the probability density distribution of the molecular wave packet. The influence of intra-molecular vibrational energy redistribution appears to be properly taken into account by the new model on the whole range of coupling strengths. It demontrates that most of the mode mixing during relaxation is due to the potential part of the Hamiltonian and not to the coupling among relaxation operators. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4892376] I. INTRODUCTION

A variety of phenomena of major importance in physics and chemistry involve the interaction between comparatively small molecules and their environment. An important example is that of chemical reactions on solid substrates, being a key step in technological processes such as catalysis, chromatography or materials processing, etc.1–10 Unravelling the atomistic mechanisms governing elementary reactions such as adsorption, bond formation and breaking, diffusion, and desorption may lead to substantial improvements in controlling the selectivity and enhancing the efficiency of the associated practical applications. A particularly interesting aspect is the role played in these processes by inter-mode coupling,11–18 which can either be mediated by the molecular potential or by weak interaction with the environmental degrees-of-freedom. In the following, we will keep the traditional use of the term “intra-molecular vibrational energy redistribution” (IVR) to refer only to the former type of coupling. Although much mechanistic information can be obtained from classical and semi-classical simulations, it is very ina) Email: [email protected]

0021-9606/2014/141(7)/074703/15/$30.00

structive to resort to a quantum mechanical description of atomic motions to gain a deeper insight into dynamical processes at surfaces. Except for a few idealized textbook examples, the solution of the quantum equations of motion describing molecules impinging on or attached to a surface must be addressed numerically.19, 20 Nevertheless, because of the unfavourable scaling of memory and computing time requirements with the increase of the dimensionality of the system, rigorous computer simulations are restricted at present to problems involving a relatively small number of dynamical variables. Accordingly, developing and testing new numerical schemes to undertake the integration of the quantum equations of motion for many-body systems continues to be crucial for achieving a realistic modelling of quantum phenomena on the nanoscale. This aspect is especially relevant in the case of molecules adsorbed on surfaces, since experiments show that the dynamics exhibits a marked correlation between the internal degrees of freedom (i.e., among the motion along and perpendicular to the surface and the vibrational and rotational motions of the adsorbate).11–18 The vibrational relaxation of a molecule in the vicinity of a non-rigid surface is mediated by the excitation of phonons21, 22 and electron-hole pairs in the solid, the second

141, 074703-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-2

L. Uranga-Piña and J. C. Tremblay

mechanism being dominant for metallic substrates.23–30 Hence, molecules scattering from a surface or sticking on it may be classified as open or dissipative systems, i.e., systems interacting with its environment. State-of-the-art quantum dynamical methods are not well suited to handle condensed phase problems because of the exponential scaling of the solution with increasing system size, but the separation of the “system” from the “bath” contributes to alleviate this difficulty. The language of open system density matrix theory provides an appealing theoretical framework to describe the reduced quantum dynamics of these phenomena, enabling to incorporate the dynamical state of the surroundings implicitly. Tracing out the environmental degrees of freedom, the dynamics along the remaining coordinates is described by the reduced density matrix (RDM), which obeys the dissipative Liouville-von Neumann equation. To map the dynamics on the reduced system, the specific dissipative model must account for the energy flow between the molecule and the environment, e.g., the creation and annihilation of phonons and electron-hole pairs in the substrate during reactions at surfaces. The phonon and the electronhole pair relaxation channels usually proceed on picosecond and femtosecond time scales, respectively.31–33 Indeed, at least for metallic hosts, the irreversible dissipative dynamics is expected to be dominated by the coupling of the molecular motion with the much faster electronic (de-)excitations of the substrate. Since the characteristic correlation time of the bath modes are frequently smaller than those of the explicitly treated degrees of freedom, short-term memory effects of the surface can be neglected and the Markov approximation holds. The resulting dissipation operator is local in time and its particular form can be derived in a variety of ways, ranging from microscopic approximate modelling of the systembath coupling to imposing physical constraints on the representability of the RDM.34–37 On the contrary, if the coupling of the excited adsorbate motion to the slow phonon modes of the substrate is significant, the bath memory time scale is not much shorter than the characteristic time scale of the system. In this case, the dissipative effects arising from the vibrational excitations in the surface cannot be disregarded. Accounting for this delayed vibrational relaxation gives rise to memory terms in the Liouville-von Neumann equation, derived from the time-correlation function of the displacements of surface atoms.38, 39 Although more sophisticated theoretical treatments have been developed,38–41 explicitly incorporating non-Markovian effects on the dissipative quantum dynamics, the influence of memory terms on the time evolution of the reduced density matrix is routinely neglected in applications of the theory of open quantum systems. Quantum systems interacting with a strong, short electromagnetic field, constitute an important example of dynamical systems lying outside the domain of applicability of the Markov approximation. The Lindblad formalism, which maps the time evolution of the full system to a Markovian semi-positive reduced dynamics, provides an analytical form for the dissipation super-operator,34 and will be used throughout this work. The strictly positive time evolution of the RDM allows for a probabilistic interpretation of its diagonal elements at all

J. Chem. Phys. 141, 074703 (2014)

times. Although the specific shape and the effect of the Lindblad operator on the dynamics is known, the rates associated with each dissipative channel must be computed from first principles. For the specific case of gas phase molecules on metallic surfaces, the theoretical framework for describing nonadiabatic interactions in the system-bath, weak coupling regime is well established.7, 42, 43 The standard procedure consists in computing the non-adiabatic couplings between the different vibronic states using first order time-dependent perturbation theory.42 The non-adiabatic coupling between electronic states mediated by the molecular vibrations arises from the nuclear kinetic energy operator. Hereafter we will consider that only first order off-diagonal vibronic couplings are relevant, therefore neglecting higher order non-adiabatic couplings. By doing so, well-established expressions for the anharmonic vibrational relaxation rates can be obtained.42, 44–47 It is important to note that the effects of IVR, i.e., of intermode coupling induced by the molecular potential, are here explicitly included in the definition of both the relaxation operators and the associate rates. Although considerable progress has been made in wave packets propagation techniques as applied to multidimensional dynamical problems,48, 49 simulating the time-evolution of the reduced density matrix according to the Liouville-von Neumann equation poses much bigger challenges in terms of the associated computational cost, as compared to the numerical solution of the Schrödinger equation.19, 20 As a rule, the feasibility of both, time-independent and time-dependent quantum mechanical calculations, relies on a judicious choice of the basis functions. Consequently, a careful evaluation of the relative advantages and shortcomings of each alternative expansion is needed in order to assess their affordability, especially for quantum systems undergoing dissipation. In the case of chemical reactions at surfaces, the reduced density matrix must often be capable of accounting for unbound motion of adsorbates, thus grid bases are the conventional choice. Nevertheless, when studying the relaxation dynamic processes of strongly coupled systems, usually confined in space, the spectral representation is more suitable, as in this case both the dissipative operators and rates are well defined. In this work, we investigate the accuracy of an approximate representation of the vibrational relaxation dynamics, which follows from three central hypotheses: (i) the system-bath coupling is weak and only first-order perturbative rates are relevant, (ii) a contracted grid basis generated via a potentialoptimized discrete variable representation (PO-DVR) can capture most of the effects of IVR in effective onedimensional grid bases, and (iii) the Lindblad operators representing each dissipative channel can be represented as a tensor product of onedimensional operators spanning the factorizable space. These assumptions result in a marked reduction of the grid size and of the scaling of the computational cost with respect to the increase of the dimensionality of the system. The physical correctness of this treatment, as compared to the reduced density matrix propagation in the basis of its

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-3

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

eigenstates, is illustrated for a two-dimensional system of two-coupled Morse oscillators with realistic parameters modelling two interacting CO molecules oriented perpendicular to a Cu(100)-surface. Two-dimensional models historically played a remarkable role in the development of the theory of molecule surface reactions50–52 and they can be regarded as a minimal representation of various phenomena such as inelastic scattering of a molecule from a substrate or the adsorption of gas phase atoms at non-rigid surfaces (i.e., the simplest model comprising the essential features of those dynamical processes). Therefore, they appear as a natural starting point to test novel theoretical approaches for which an extension to larger dimensionality problems is envisaged. In Sec. II, the reduced density matrix formalism and the perturbative separable scheme employed to model the adsorbates coupling to the internal modes of the substrate are presented. Sec. III introduces the model system under study as well as the computational details regarding the time integration. In Sec. IV, the numerical procedure used to generate the PO-DVR of the vibrational motion is introduced. The relaxation dynamics of the system is analysed in Sec. V, in terms of the time evolving state populations and the density distribution of the molecular wave packet. The conclusions are briefly summarized in Sec. VI.

B. Choice of dissipative operators

There are several choices for selecting the explicit form of the bath operators Cˆ mn (t) depending on the physical effect to be described. For instance, in the case of pure dephasing, i.e., when the bath cause decoherence in the system, but no energy is exchanged, the Lindblad operators are diagonal in the basis of the system eigenstates: 1/2 Cˆ mn = δmn n→m |mm|.

represent environment-induced The quantities transition rates between states n and m, and they are computed perturbatively as described bellow. Since dissipative operators Cˆ mn are diagonal, the populations ρ mm (t) remain unchanged as the system evolves, thus conserving the energy during the evolution. Moreover, the coherences given by the off-diagonal terms of the density matrix, ρ nm (t), decay exponentially. This scenario describes the case where elastic collisions between system and bath particles take place. The other limiting case occurs when the system relaxes via inelastic collisions with the environment, while also undergoing dephasing. In this case the Lindblad operators are strictly non-diagonal, i.e., 1/2 Cˆ mn = n→m |mn|, m = n,

(4)

Cˆ mn = 0, m = n.

(5)

II. QUANTUM DISSIPATIVE DYNAMICS 1/2 n→m

stands for the rate of population transfer from Again, state n to state m.

A. General Liouville-von Neumann equation

Open quantum systems exhibit irreversible time evolution due to their interaction with the environment. Therefore, a reduced description of the system dynamics is required. Assuming that the total evolution of the system together with its environment is unitary, then the dynamics of the reduced system can be mapped to follow the evolution of a completely positive semigroup.34, 53 Thus, the time evolution of the re˙ˆ of a Markovian system is governed duced density matrix ρ(t) by the Liouville-von Neumann equation, given here in the Schrödinger representation: i ρ˙ˆ = − [Hˆ s , ρ] ˆ + LD ρ, ˆ ¯

(1)

where Hˆ s is the system Hamiltonian and LD is the Lindbladian superoperator that accounts for dissipation effects. The explicit form of this superoperator is LD =

1 mn

2

(3)

1/2 n→m

† † ([Cˆ mn (t), ρ(t) ˆ Cˆ mn (t)] + [Cˆ mn (t)ρ(t), ˆ Cˆ mn (t)]).

(2) Here, Cˆ mn (t) are the coupling operators accounting for the interaction between the system and the bath, in particular representing the (m → n) dissipative channel. The bath operators Cˆ mn can be derived microscopically starting from the systembath Hamiltonian, taking the strict secular limit of the secondorder Born-Markov approximation to the time-evolution of the reduced density matrix.36, 54 In general, the operators Cˆ mn will describe the “jumps” between eigenstates of Hˆ s induced by the interaction with the bath.

C. Transition rates

As seen previously, for energy relaxation the spectral basis, Hˆ s |n = En |n,

(6)

appears as the natural choice for expressing the Lindblad operators and, hence, for the equations of motion (EOM). In the interaction picture and in the spectral basis representation, this EOM read  (I ) (t) = δmn n →m ρn(I n)  (t) ρ˙mn n



  n

m→n

+ n→n 2

 (I ) ρmn (t).

(7)

The transition rates  n→m mix the coupled multidimensional states |n and |m via the system-bath coupling Hamiltonian. As it has been demonstrated in Ref. 42, proper inclusion of the inter-mode coupling induced by the molecular potential can lead to drastically different microscopic relaxation scenarios. In the weak-coupling limit and factorising the vibronic states, the transition rates can be evaluated using Fermi’s golden rule. At zero temperature it reads 2π  |n|i|Wˆ |f |m|2 δ(Ef − Ei + ¯ωmn ), (8) n→m = ¯ if where Wˆ is the interacting system-bath Hamiltonian. Ei and Ef are all possible initial and final state energies of the bath

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-4

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

modes whereas the bath energy uptake (loss), (Ef − Ei ), is transferred from the system into the surroundings allowing the former to relax. Within the context of an adsorbate approaching a metallic surface, the mechanism of energy transfer between the vibrational degrees of freedom of the molecule and electronic degrees of freedom of the surface is well known.7, 42–44, 46 Essentially, the electron-hole pairs creation in the metal is associated with the transfer of population between two vibrational states of the molecule. In this particular case, the system-bath coupling is nothing else but the kinetic energy operator of the adsorbate nuclei. The transition rates are then expressed as 2 2   2π     − ¯ n|i|∇q2 |f |m n→m =   ¯ if 2Mq q  × δ(Ef − Ei + ¯ωmn ) =

 π ¯3  |(i|∇q2 |f n|m 2 2M q q if + 2i|∇q |f n|∇q |m)|2 δ(Ef − Ei + ¯ωmn ).

(9)

Here, Mq is the mass corresponding to the motion along the qth degree of freedom. Note that, in the spirit of first order perturbation theory, the contribution of each mode q to the transition rate is treated additively as a separate channel. The first term of Eq. (9) involves second derivatives of the electronic wave functions and is usually much smaller than the second term involving only first derivatives. Indeed, the former can be rewritten as second-order contributions and they will be neglected in the present treatment. The approximate form of the transition rates thus reads   (q) |n|∇q |m|2 , (10) n→m =

number of dissipation channels to be considered, solving the secular equation (6) may require a great effort from the computational standpoint. This constitutes the main drawback of this particular choice of Lindblad operators when considering multiple channels and several degrees of freedom. The main contribution of this article is to propose a more suitable form for the dissipative operators, in order to improve the former choice and, as a consequence, to make numerically converged calculations feasible for larger systems. D. Alternative representation of Lindblad operators

Given the particular form of the system-bath Hamiltonian as a sum of kinetic energy operators corresponding to each nucleus, and therefore the corresponding structure of the transition rates between two eigenstates |n and |m, it becomes natural to think about exploring an alternative representation of the Lindblad operators using a more suitable basis set. We propose to construct the multidimensional dissipative operators as a tensor product of one-dimensional operators corresponding to each degree of freedom. The question that arises right away is whether this representation is sufficient to correctly describe the intra-molecular couplings, that were explicitly included in the former representation within the multimode states |n. Hereinafter, we will answer this question by solving the resulting equation of motion for a model system. The alternative choice of dissipative operators results in transition rates which are naturally decomposed into a sum of corresponding one-dimensional rates associated with each degree of freedom. The Lindblad operators in the tensor product basis set read 1/2 Cˆ m n = n0 →m0 |m0 n0 |,

where

q

where the electronic contributions are given by Refs. 42 and 44 2π ¯3  i|∇q |f 2 δ(Ei − εF )δ(Ef − εF ). (11)  (q)  Mq2 if The definition of the scaling constant is obtained in the quasistatic limit, where ¯ωmn → 0 and the initial and final electronic states are about the Fermi energy εF . This approximation does not affect the rates at larger transition frequencies, since the integrals on the right-hand side of Eq. (10) vanish as the energy difference increases. It is important to stress that, even after obtaining a simplified expression for the transition rates (Eq. (10)), constructing the Lindblad operators in the representation of the system eigenstates, as in Eq. (4), requires the knowledge of the multidimensional nuclear wave-functions, |n. In practice, the requirement to perform the computation of the  n→m effectively limits the amount of degrees of freedom that can be treated, especially if a great number of dissipative channels have to be considered. Indeed, it restricts the actual possibilities of numerically solving the equation of motion (7), owing to the fact that obtaining the solutions of Eq. (6) is a prerequisite for the evaluation of the transition rates. Depending on the

(12)

0 0

|m0  = |mq  ⊗ |mq  · · · ⊗ |mq , 1

2

(13)

N

N being the total number of degrees of freedom of the system. Furthermore, the new transition rates n →m can be split into 0 0 N one-dimensional rates, each one representing the singlemode dissipative channels. Within this approach, the transition rates become decoupled: m →n = 0

N 

0

Iq ⊗ Iq ⊗ . . . ˜ m 1

i

2

q

i

→nq

⊗ . . . Iq . N

i

(14)

Henceforth, the advantages of this assumption to tackle the quantum reaction dynamics of systems of sequentially increasing dimensionality are striking: the numerical effort associated with the computation of the quantities m →n scales 0 0 linearly with the system dimension. If the system Hamiltonian is conveniently decomposed into a sum of one-dimensional Hamiltonians plus an intra-molecular coupling term:  Hˆ = Hˆ 0 + Vˆ , Hˆ 0 = Hˆ q , (15) i

i

then the states |m0  can be chosen as the eigenfunctions of the problem Hˆ 0 |m0  = Em |m0 . In that case, the 0 Liouville-von Neumann equation in the interaction picture

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-5

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

becomes

III. MODEL SYSTEM

i (I ) ρ˙ˆ (I ) (t) = − [Vˆ , ρ (I ) (t)] ¯   mq →nq   † i i Aˆ m n (t), ρˆ (I ) (t)Aˆ mq nq (t) + q q i i 2 i i i m n

This section is devoted to testing the new approach in a two-dimensional system interacting with a bath. The Hamiltonian of the model system was chosen to resemble the realistic interactions of two coupled CO molecules approaching a Cu(100) surface. The parameters of the interaction potentials between the CO molecules and the copper substrate are adapted from Refs. 56 and 57. Previous calculations carried out for the CO/Cu system reference 39, including both instantaneous and delayed dissipation, have shown that the main effect of the latter is the appearance of quantum coherences in the time evolution of the state populations. Nevertheless, the overall population dynamics remains similar to that obtained if only instantaneous dissipation is considered. Delayed dissipation is expected to play a minor role in the present case since only normal vibrations of the CO centre of mass relative to the surface are considered here, and the surface vibrations couple preferentially to the frustrated translational mode. Therefore, in this paper, we consider Lindblad instantaneous rates only.

q



+ Aˆ m

q

where Aˆ m

q

i

nq

i

q

i

nq

i

 † (t) ρˆ (I ) (t), Aˆ mq nq (t) , i

i

(16)

i

= |mq nq |. The interpretation of this equai

i

i

tion of motion becomes clear: while the alternative Lindblad operators accounts for the dissipation channels of onedimensional modes coupled to the bath, the first term of the equation containing Vˆ , properly incorporates the intramolecular vibrational energy redistribution between the system degrees of freedom. That is, by choosing a representation for the dissipative operators of uncoupled uni-dimensional channels linked with the single-mode Hamiltonians, we are thus transferring the coupling between the system degrees of freedom into the commutator between the density matrix and Vˆ . When Vˆ = 0, both representations coincide. Care must be taken in choosing the basis set for the representation of the operators. Since the dynamics of the system is not unitary, it evolves towards a pointer state that depends on the chosen basis.36 In the Lindblad form, proper unitary transformation of all operators will yield the same microscopic relaxation scenario at the cost of losing the diagonal form of the dissipative super-operator.36, 55 Using a tensor product basis of one-dimensional pseudo-spectral functions, |m0 , the equations of motion take the form i  (I ) (I ) Vm(I )n ρn(I n) (t) − ρm  (t)Vn n n 0 0 0 0 ¯ n  + δm n n →m ρn(I n)  (t)

(I ) ρ˙m n (t) = − 0 0

0 0



0

n

  m →n + n 0

n

2

0

→n

 (I ) ρm n (t). 0 0

(17)

The matrix exhibits lots of structure, which can be exploited to render the time integration of Eq. (17) very efficient while simultaneously limiting the memory storage of the Liouville super-operator. Despite the apparently more complicated structure of Eq. (17) as compared to Eq. (7), the former is, indeed, more suitable when the number of degrees of freedom increases, as solving Eq. (7) becomes infeasible. Once again, the advantage of Eq. (17) over Eq. (7) lies in the possibility of computing the high dimensional transition rates as a tensor product of one-dimensional rates, via the one-dimensional eigenstates, |mq . The costly diagonalization of the highi dimensional Hamiltonian, Hˆ , to generate the full spectral basis and to compute the high-dimensional transition rates, is no longer feasible when the number of degrees of freedom increases. The grid-based representation appears thus more computationally suitable for high-dimensional quantum dynamics, especially if the stochastic unravelling of the RDM is preferred to the full density matrix propagation technique.

A. Hamiltonian

In Fig. 1 a cartoon of the problem we use as a toy system is depicted. The system Hamiltonian for two rigid CO molecules oriented perpendicular to the copper surface is given by Hˆ = Tˆq + Tˆq + Vˆ1 (q1 ) + Vˆ2 (q2 ) + V cpl (q1 , q2 ), 1

2

(18)

where q1 and q2 are the respective distances from each of the CO-centre of masses to the surface. Tˆq and Tˆq are the one1 2 dimensional kinetic energy operators while Vˆ1 , Vˆ2 , and V cpl designate, respectively, the potential energy operators corresponding to the motion along the q1 and q2 coordinates and the coupling between these two vibrational modes. Here we adopt the dissipative model introduced in Ref. 57, where stateresolved anharmonic transition rates were calculated perturbatively, as in Eq. (10). The non-adiabatic couplings  m→n , between the adsorbates degrees of freedom and the electrons at the surface, drive the population transfer from state |m to state |n. Upon relaxation, dephasing destroys coherence in the system.

FIG. 1. Scheme of the system: two CO molecules approaching a Cu(100) surface. q1 and q2 are the distances from each CO centre of mass to the surface, respectively. d corresponds to the distance between the centre of mass of each molecule.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-6

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

B. Potential energy surface

The Hamiltonian Hˆ (Eq. (18)) describing the vibrations of the target two-dimensional system comprises the onedimensional potential energy operators Vˆ1 (q1 ) and Vˆ2 (q2 ), which are taken as Morse functions with equilibrium diseq tance qi = 4.805 bohrs, well depth Dq = 0.02128 a.u., and i anharmonicity parameter αq = 1.588 bohr−1 , (i = 1, 2). i The parameter values are chosen so that the resulting Morse potential resembles the interaction between a CO molecule and a Cu(100) substrate extracted from Ref. 56. On the other hand, the non-separable term V cpl (q1 , q2 ) consists in a screened bi-linear coupling of the form: ⎡ ⎤ Mq Mq ωq ωq 1 2 1 2 ⎦ (q1 − q1eq )F1 (q1 ) V cpl (q1 , q2 ) = a ⎣ 2 eq

· (q2 − q2 )Fq (q2 ),

(19)

2

where ωq designate the oscillation frequency (within the hari monic approximation) along each coordinate. This magnitude of the prefactor takes the value of 0.00086 a.u. The screening functions F1 , F2 guarantee the linear behavior of the coupling in the neighbourhood of the equilibrium position along each degree of freedom and they asymptotically vanish for large elongations. In this work, the applicability of the numerical scheme described above, to account for the dissipation dynamics of two coupled identical Morse oscillators, is explored for several coupling strengths. The dimensionless parameter a, which controls the relative intensity of the interaction term V cpl (q1 , q2 ) with respect to the uncoupled Hamiltonian, is varied within the range from 0.0 to 1.0. It can be understood as the distance-dependence (d) of the interaction between the two CO molecules. The value a = 1 refers to a situation where the coupling and the one-dimensional potentials Vˆq are of the same order of magnitude. As a consei quence of the symmetry of the problem, the screening functions are also taken to be equal:  q 8 ⎡  q 8 ⎤ Fi (qi ) = e



i eq qi +2/α

q

i

· ⎣1 − e



i eq qi −1/α

q

i

⎦.

(20)

the direct diagonalization of the two-dimensional Hamiltonian Hˆ is necessary. To this aim, Hˆ is represented using a tensor product basis built up by the one-dimensional potential optimized-discrete variable representation (PO-DVR) basis sets corresponding to each degree of freedom, |α|β. The two-dimensional Hamiltonian eigenfunctions can be then written as  Cαβ,n |α|β. (21) |n = αβ

In this basis, the matrix elements of the Hamiltonian exhibit also a favourable sparse structure: (1) (2) Hαβ,α β  = Tα,α  δββ  + Tβ,β  δαα  + (V1 (q1α ) + V2 (q2β )

+ V cpl (q1α , q2β ))δαα δββ  .

In the PO-DVR basis, the kinetic energy operators are onedimensional matrices and the potential is diagonal and simply evaluated at the DVR points, α, β. The diagonalization of the Hamiltonian Hˆ in that basis yields the expansion coefficients Cαβ, n . In a second approach, we split the Hamiltonian Hˆ as follows: Hˆ = Hˆ 0 + Vˆ .

C. Basis representation

The calculation of the derivative in Eq. (1) implies the evaluation of the action of both the Liouvillian and the Lindblad superoperators on the density operator at every time step, which can be expressed more conveniently in the Dirac representation. Accordingly, two different routes were followed to describe the vibrational relaxation of the system: first, by solving the Liouville-von Neumann equation (1) in the full spectral basis representation (7), and second, by using the direct product basis representation corresponding to the effective one-dimensional Hamiltonians leading to Eq. (16) and then projecting on the same basis (Eq. (17)). In the first case,

(23)

The contribution,

  Hˆ 0 = Hˆ qeff + Tˆq + Vˆqeff , + Hˆ qeff = Tˆq + Vˆqeff 1 2 1 2 1

2

(24)

is a simple sum of effective one-dimensional Hamiltonians, and Hˆ qeff , generated by the PO-DVR procedure which is Hˆ qeff 1 2 described bellow. The advantage of using this basis set lies in the reduction, to a large extent, of the intra-molecular coupling originally present in Hˆ . The residual potential coupling, Vˆ , is given by   + Vˆq − Vˆqeff . Vˆ = Hˆ − Hˆ 0 = Vˆ cpl + Vˆq − Vˆqeff 1 2 1

2

(25) This partition enables the use of the effective Hamiltonians eigenbasis: Hˆ qeff |mq  = q |mq , (i = 1, 2), i i

The analysis of the effect of the potential asymmetry on the numerical performance of the proposed method is deferred for further discussion in a forthcoming contribution.

(22)

i

i

(26)

to perform the transformation to the interaction picture. The solution of Eq. (1) requires the specification of an initial state, generally mimicking certain experimental conditions. Here, the initial density matrix corresponds to a pure vibrational state along each degree of freedom, which is denoted as |2q ⊗ |0q in product representation, or as (2, 0) in 1 2 pseudo-spectral form. D. Numerical details

A fourth-order Adams-Moulton predictor-corrector method, initiated by the fourth-order Runge-Kutta integrator, was used to evolve in time the matrix elements of the reduced density matrix according to the dissipative Liouville-von Neumann equations (7) and (17). Convergence tests were carried out for the different parameters influencing the numerical propagation, e.g., the time step and the number

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-7

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

of basis functions in the grid representation introduced in Sec. III C, and a step size of 5 a.u. was chosen. The total propagation time was set to 24 ps. For comparison, the lifetime of the first excited state along the mode q1 was calculated from first principles in the local harmonic approximation. It was found to be 3 ps, in good agreement with experimental estimates. At each time step, the projection operators Pˆn = |nq  q 1 1 nq | ⊗ Iq and Pˆn = Iq ⊗ |nq nq | are used, respectively, to compute the state populations pn (t) and pn (t) according q q 1 2 to

calculations is constructed as the tensor product of the eigenfunctions |mq  for each degree of freedom qi . i The definition of the effective Hamiltonians in (26) is also not unique. However, a clever partitioning of the total Hamiltonian may help to shift large parts of the coupling into the effective operators Hˆ qeff . In doing so, the commutator appearing i in Eq. (16) containing the coupling gets reduced, ultimately becoming negligible and yielding a purely dissipative dynamics in the interaction representation. Following the original implementation by Echave and Clary,58 we define the effective contributions entering the terms in the right-hand side of Eq. (23) as

pn (t) = Pˆn  = Tr{Pˆn ρ}, ˆ

Vˆ1eff (q1 ) = Vˆ1 (q1 ) + V¯ (q1 ), and

q

q

1

2

2

q

1

2

q

1

(27)

1

Vˆ2eff (q2 ) = Vˆ2 (q2 ) + V˜ (q2 ),

with similar formula holding for the populations pn (t). Ad2

ditional observables were calculated to aid the interpretation of the physical process under study, for instance, the purities, (28)

as well as the time-averaged mean squares of the onedimensional-populations and -purities differences,   2 1 T  f ull f act pnq − pnq pi = dt, (29) i i T 0 and

 γ =

1 T



T

(γ f ull − γ f act )2 dt,

(30)

0

where T is the total propagation time. The superscripts fact and full stand for the factorization and full-dimensional procedures used to compute the transition rates. IV. POTENTIAL-OPTIMIZED DISCRETE VARIABLE REPRESENTATION (PO-DVR)

As it is standard in discrete variable representations, the matrix elements of differential operators such as the kinetic energy or the transition rates are computed exactly, while those operators which are local in the position operator representation (e.g., the potential energy contributions) are evaluated approximately, within the accuracy of the Gaussian quadrature. The PO-DVR is based on the use of an underlying dense grid to solve the time-independent Schrödinger equations:  |mq  = q |mq , (i = 1, 2), (31) |mq  = Tˆq + Vˆqeff Hˆ qeff i i i

i

i

i

i

which provides numerically exact one-dimensional eigenfunctions for the low-energy states. Let us recall that we are neglecting the inter-mode coupling in the calculation of the one-dimensional eigenstates. The set of eigenfunctions |mq  i is filtered by imposing certain threshold max on the eigenvalues q . This filter is not unique but allows for a good descripi tion of low-energy states, which is relevant to the dynamics reported in this work. In a second stage, the resulting contracted basis sets are used to represent the one-dimensional effective Hamiltonians, Hqeff , which are subsequently diagi onalized. The representation actually used in the dynamical

V¯ (q1 ) = min(V cpl (q1 , q2 )),

(32)

V˜ (q2 ) = min(V cpl (q1 , q2 ) − V¯ (q1 )).

(33)

q2

q1

Here V¯ (q1 ) is the coupling between the Morse oscillators obtained by minimizing, for each value of the coordinate q1 , the superposition Vq (q2 ) + V cpl (q1 , q2 ) with respect to co2 ordinate q2 . The same procedure is used to generate V˜ (q2 ), where the minimization is now performed with respect to coordinate q1 . In this case, the second term, k1 = min(V¯ (q1 )), is q1

just a constant. Within this partitioning, the expression for the residual potential term Vˆ as a function of V cpl (q1 , q2 ) takes the following form: Vˆ (q1 , q2 ) = V cpl (q1 , q2 ) − V¯ (q1 ) − V˜ (q2 ).

(34)

In this paper, this quantity is varied to evaluate the equivalence between the factorizable and the fully coupled ansatz, allowing for a systematic, unbiased investigation of the effect of the IVR over a wide range of intermode coupling strengths. In Fig. 2, the ratio between Frobenius norms corresponding to matrix elements representing the operator Vˆ and those of the potential energy operator Vˆ1 (q1 ) + Vˆ2 (q2 ) 0.01

1e-07

H

2

V

2 ΔVij

2

f = Σ ΔVij / Σ Hij

0.008

f =Σ

V

γ (R) = Tr (ρˆ (R) ρˆ (R) ), R = {f act, f ull};

where the auxiliary functions V¯ (q1 ) and V˜ (q2 ) are given by the expressions:

Frobenius norm (f )

q

/ Σ

8e-08

2 Vij

H

2

0.006

6e-08

0.004

4e-08

0.002

2e-08

0

0

0.2

0.4 0.6 coupling strength (a)

Frobenius norm (f )

1

0.8

1

0

FIG. 2. Black points: ratio between Frobenius norms of matrices Vij and Vij for different inter-mode coupling strengths. Red line: ratio between Frobenius norms of matrices Vij and Hij for different values of a.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-8

L. Uranga-Piña and J. C. Tremblay

+ V cpl (q1 , q2 ) are shown. The mathematical expression for this ratio is  2 i,j |Vij | V f =  . (35) 2 i,j |Vij | It can be observed, that the residual coupling Vˆ represents only a small fraction of the total interaction potential. The same behavior is observed for the ratio between Frobenius norms of matrix elements corresponding to Vˆ and the total energy of the system, respectively:  2 i,j |Vij | H f =  . (36) 2 i,j |Hij | In both cases these magnitudes have values below the 1%, thus indicating that the residual potential operator Vˆ represents indeed a relatively small perturbation. Fig. 3 shows a comparison between the coupling V cpl (q1 , q2 ) and the residual potential term (Vˆ (q1 , q2 ) − k1 ), obtained via the transformation procedure outlined above. The constant term, k1 , contained in Vˆ has been disregarded. As it can be seen, the procedure employed tends to flatten the attractive part of the coupling potential towards zero, while increasing the repulsive part. This is a result of the filtering choice, which transfers an important part of the . coupling to the effective one-dimensional Hamiltonians, Hqeff i

J. Chem. Phys. 141, 074703 (2014)

This behavior tends to be more pronounced as the coupling strength increases, and it can be observed more clearly if we look at the region in which the dynamics takes place (inside the dashed square). This area ranges approximately from 4.6 Å to 5.2 Å along both axes and, as it can be seen, the flattening effect of the coupling potential is more marked for the largest coupling constant a = 1.0 (bottom panels). Upon construction, the Hamiltonian Hˆ 0 contains a significant part of the correlation between the two degrees of freedom in the region relevant for the dissipative dynamics and an efficient and most likely accurate first approximation would be to disregard the potential coupling term in the Liouvillevon Neumann equation of motion. Still, our goal is to assess the influence, on the coupled dynamics, of the choice of the ansatz (factorizable or full) introduced in order to evaluate the transition rates. Therefore, the results presented in Sec. V were obtained by considering both diffusion and dissipation terms in Eqs. (7) and (17).

V. RESULTS AND DISCUSSION

As mentioned above, the central question to be addressed is whether expressing the Lindblad operators as a tensor product of one-dimensional dissipative operators, and thus transferring the coupling between the degrees of freedom from the two-dimensional rate constants towards the commutator part of Eq. (16), is sufficient to recover the same system dynamics derived from solving the full two-dimensional Liouville-von Neumann Eq. (7) with two-dimensional Lindblad operators. To assess the influence of the product ansatz, in this section we will present the comparison between different observables calculated within these two approaches that we shall call hereafter the factorizable and the full methods, respectively, for the specific system described in Sec. III. First, we will examine the time evolution of the onedimensional state populations corresponding to each degree of freedom, pn (t) and pn (t), as introduced in Eq. (27). q

q

1

2

Fig. 4 shows the behavior of these quantities during the first 24 ps of the dynamics. Taking into account the initial conditions and that exclusively downward transitions are allowed, at T = 0, according to the principle of detailed balance, only those states fulfilling nq , nq ≤ 2 are considered. The various 1 2 panels, from top to bottom, correspond to different coupling constants ranging from the smallest a = 0.0 to the largest one a = 1.0. The left panels correspond to the population pn (t) q

1

of states nq ∈ {0, 1, 2}, while the corresponding state popu1 lations pn (t) are displayed on the right. q

FIG. 3. Comparison between the coupling potential V cpl (q1 , q2 ) (left panels) and the residual potential V (q1 , q2 ) (right panels) after applying the PO-DVR procedure, for different inter-mode coupling strengths. The region where the dynamics takes place lies between 4.6 Å and 5.2 Å approximately, along both axes. The resulting coupling potential V tends to flatten around this region compared with V , even for the largest inter-mode coupling a = 1.0.

2

For the case of vanishing coupling strength (a = 0.0), the population dynamics serves as a test of the numerical stability of the algorithm, as both approaches have to yield same results for this limiting case. Since there is no population transfer between the two degrees of freedom, the probabilities pn (t) (nq = 0) remain equal to zero at all times and q

2

2

the problem is effectively one-dimensional. As for pn (t), the q

values pn (0) = δ2 n q

1

q

1

follow from the initial conditions. It 1

can be seen that the state nq = 2 undergoes an exponential 1

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

L. Uranga-Piña and J. C. Tremblay

pn

a = 0.0

pn

a = 0.1

q1

q1

pn

J. Chem. Phys. 141, 074703 (2014)

n = 0 (full) n = 0 (fact) n = 1 (full) n = 1 (fact) n = 2 (full) n = 2 (fact)

a = 0.3

q1

1 a = 0.0

pn

a = 0.1

pn

a = 0.3

a = 0.1

q2

q2

pn = 0 q pn = 0 q pn = 1 q pn = 1 q pn = 2 q pn = 2

Population

074703-9

1

2

1

2

pn

1

q2

q2

pn

a = 0.5

q1

a = 0.5

pn

q2

0

0

5

10

15

20

t [ps] pn

a = 1.0

q1

0

5

10

15 t [ps]

20

a = 1.0

0

5

10

15

pn

q2

i

q

20

t [ps]

q

q

1

and pn , corresponding to the inter-mode coupling strength a = 0.1. Notice

FIG. 4. Comparison between the effective one-dimensional state populations pn calculated by both, the full and factorizable methods. The pn correq

FIG. 5. Time-evolution of the effective one-dimensional populations pn

i

spond to the population of each degree of freedom qi (i = 1, 2) and each effective single-particle state nq = 1, 2, 3, for different coupling strengths a. i The results obtained applying the factorizable method converge quite well to those obtained from the full approach even for the largest inter-molecular coupling. The oscillation frequency increases with the coupling strengths (see explanation within the text).

decay. The population of both the ground and the first excited vibrational states increase during the first 5 ps, the growth of the latter taking place at a faster rate. From 5 ps onwards, the population of the state nq = 1 is also progressively depleted 1 as a consequence of the population transfer to the ground state. The calculations carried via the bi-dimensional and the factorizable representations of the dissipation operators yield identical results, as expected. The overall relaxation of the system takes place within the first 24 ps approximately. This close correspondence between the two simulation schemes persists upon inclusion of non-vanishing couplings. For a > 0.0, the relaxation stops to be restricted to a single vibrational mode but the predictions based on the factorization method (Eq. (17)) stay in excellent agreement with the dynamics derived within the full-dimensionality approach. For a = 0.1, it can be noticed that the population dynamics roughly follows the same trend as the uncoupled case, namely the decrease of the population of the initial state nq = 2 and 1 the gradual increase of the probabilities to find the system in states nq = 1 and nq = 0. However, for this coupling 1 1 strength, a marked oscillatory pattern gets superimposed on the average behavior of the state populations. While the average trend of the one-dimensional state populations is a fingerprint of the effect of dissipation, the oscillations arise from the inter-mode coupled dynamics, more specifically from the IVR. Thus, the deviations from the former becomes increasingly more pronounced as the coupling gets larger. The oscillatory pattern exhibited by the effective onedimensional populations can be understood if we consider that after tracing out one of the degrees of freedom, the intermode coupling can be regarded as an external oscillatory field

2

the overlapping between states nq = 1 and nq = 1, whereas states (nq = 2 1 2 1 and nq = 2), and (nq = 0 and nq = 0) have opposite phases. This is a fin2 1 2 gerprint of the presence of the vibrational coherence in the evolution of the density matrix.

affecting the dynamics of the remaining vibrational mode. This scenario resembles the problem of a three-level system subject to an oscillatory field with constant and equal Rabi frequencies, in which the Rabi frequencies are smaller than the detuning of the laser field.59, 60 To further understand this behavior let us focus on the evolution of the population pn corresponding to the coupling strength a = 0.1, q

1

as depicted in Fig. 5. From the figure it can be seen that populations corresponding to states nq = 0 and nq = 2 are 1 1 completely correlated and with opposed phases. The oscillatory pattern superimposed to the general dissipative trend has clearly two frequencies, whereas the one corresponding to nq = 1 has only one. The population flows initially from 1 state nq = 2 to nq = 0 passing trough state nq = 1 and re1 1 1 turning, performing further periodic cycles. The population in state nq = 1 comes from both contributions, namely in tran1 sit from nq = 2 → 1 → 0 and vice versa nq = 0 → 1 → 2, 1 1 leading to a beating frequency twice that of the fastest oscillatory component corresponding to the transition nq = 0 → 2. 1 On the other hand, Fig. 5 shows how the populations in q1 -axis compare to those of the q2 -axis. It can be clearly seen that the states nq = 1 and nq = 1 completely over1  2 lap whereas both pair of states nq = 2 and nq = 2 , and 1 2  nq = 0 and nq = 0 have opposite phases. The behavior of 1 2 the nq = 1 and nq = 1 states is a fingerprint of the presence 1 2 of vibrational coherence in the time-evolving density matrix. In other words, IVR remains a coherent process in the coupling range investigated here. By further increasing the parameter a up to the value of 0.3, an interesting behavior is observed: the curves representing the populations of the ground and the second excited state become appreciably more correlated, pointing to the more important role played by the transitions between these two states mediated by the intermediate one. Furthermore, as the coupling strength parameter is increased, the oscillations of the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

L. Uranga-Piña and J. C. Tremblay

full

fact. (t = 0.48 ps) q 2 [a.u.]

5.2

J. Chem. Phys. 141, 074703 (2014)

5.0 4.8 4.6

(t = 1.92 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

(t = 5.52 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

5.2

(t = 10.80 ps) q 2 [a.u.]

(t = 10.80 ps) q 2 [a.u.]

(t = 6.00 ps) q 2 [a.u.]

(t = 3.60 ps) q 2 [a.u.]

(t = 0.00 ps) q 2 [a.u.]

074703-10

5.0 4.8 4.6 4.6 4.8 5.0 5.2 q 1 [a.u.]

4.6 4.8 5.0 5.2 q 1 [a.u.]

(a)

5.2

full

fact.

4.6 4.8 5.0 5.2 q 1 [a.u.]

4.6 4.8 5.0 5.2 q 1 [a.u.]

5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

(b)

FIG. 6. Snapshots of the evolution of the probability distribution for certain points in time calculated with both approaches: full (left panels) and factorizable (right panels). The initial condition is set as (2, 0). (a) Zero inter-molecular coupling. The dynamics takes place only in one axis and both approaches yield identical results in this limiting case. (b) Intermediate coupling (a = 0.5). Strong IVR causes fast energy flow, leading to the inversion of population along each degree of freedom. Both approaches yield similar results.

different populations become more and more pronounced as it is expected, due to the larger Rabi frequency. Even for the strongest coupling, the population evolution predicted by the two methods agrees remarkably, thus pointing to the conclusion that the one-dimensional tensor product ansatz of the multidimensional dissipative operators provides a quite satisfactory description of multidimensional quantum relaxation processes. Figs. 6(a) and 6(b) depict the evolution of the different wave packets along both q1 and q2 axis, corresponding to distinct coupling strength parameters and to both computational approaches employed in this work (similar curves for other coupling strengths can be found in the Appendix, Figs. 9(a)– 9(d)). The representative snapshots were chosen to highlight important changes in the nodal structures along the dissipative reaction path. Fig. 6(a) shows the density evolution corresponding to certain points in time for the case of uncoupled oscillators. It can be observed that, first, the dynamics takes place only along the q1 -axis as it is expected taking into account the initial conditions. That is, the initial population is concentrated in nq = 2 and nq = 0 states and de1 2 cays smoothly towards state nq = 0 via the intermediate level 1 nq = 1 (3.60 ps) within the first 11 ps. 1 Moreover, it can be clearly noticed that both methods behave identically, as it must be for the limiting case of zero coupling. By further increasing the coupling parameter we arrive to the scenario depicted in Fig. 6(b) corresponding to a = 0.5, the intermediate coupling strength, in which the two modes are strongly influencing each other and the strong

IVR causes fast energy flow between them. The energy exchange leads to the inversion of population along each degree of freedom, compared to the initial state, i.e., when going from (2, 0) to (0, 2) (the notation (nq , nq ) was chosen to 1 2 label the coupled two-dimensional eigenstates). This particular feature is present at several points in time and it occurs more frequently as the coupling parameter increases, as it can be noticed from Fig. 4 in the spectral representation, and in Fig. 6(b) in the grid representation for the specific time instants t = 0.48 ps and t = 1.92 ps. After switching on the coupling, the two-dimensional wave packet can not be described anymore as a Hartree product of single-particle wave functions corresponding to each dynamical variable, but as a linear superposition of all possible combinations of Hartree products corresponding to different energy states of each mode. As a consequence, the wave packet turns out to exhibit the intricate shape observed, for instance, at t = 5.52 ps and t = 10.80 ps in Fig. 6(b). Additionally, it can be confirmed that for intermediate couplings the full and the factorizable approaches converge approximately to the same values. For systems subject to strong IVR, the overall dynamics consists on the periodic, coherent broadening and contraction of the probability density along the two coordinates, eventually relaxing to the same ground state wave function (0, 0). Even at the largest coupling strength, the two methods behave almost identically (see the Appendix), thereby reinforcing the overall conclusion that the factorizable method entails the same dynamics as the full for all possible coupling strengths, provided prior effort is made to include as much of the coupling in the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

Purity

074703-11

a = 0.0

a = 0.1

a = 0.3

a = 0.5

a = 0.8

a = 1.0

(full)

Purity

γ

(fact)

Purity

γ

0

5

10 15 t [ps]

20

0

5

10 15 t [ps]

20

FIG. 7. Time-evolution of the purity for each of the coupling strengths at T = 0 K. All purities start at 1 since the initial condition (2, 0) corresponds to a pure state in all cases. The initial loss of purity is associated to the initial decoherence due to the interactions between the system and the bath particles. As the system relaxes, it reaches its ground state and the purity approaches towards the unity value. The oscillations of the purity within the full approach compared with the straight behavior of the factorizable approach, can be explained by the choice of the initial condition (see text). Both approaches yield very similar results.

effective one-dimensional operators defining the pseudospectral basis. This is to say that potential-mediated IVR is the dominant effect determining the intermode mixing upon relaxation, while the coupling between relaxation operators remains only marginal. Thus, it is possible to utilise the advantageous form of the tensor product of one-dimensional dissipative operators to perform high-dimensional quantum calculations, shifting the correlation between these single-mode wave functions acting as a basis set, to the coupling part of the Hamiltonian. The evolution of the system purity for each value of the coupling strength parameter as well as for each of the approaches considered in this article (factorizable and full), is depicted in Fig. 7. Since the system is initially prepared as a coherent superposition of the zeroth-order states of Hˆ , the initial state has a purity equal to 1, independent of the basis choice. Very quickly, the purity starts to decrease from the value 1 till reaching a minimum in every case, this phase corresponding to the initial interaction between the system and the bath particles. The overall trend of the time evolution of the purity can be divided into two time steps: the initial decreasing behavior, and the later increase in magnitude towards the asymptotic limit of unity. The first contribution is due to the loss of quantum coherence of the wave packet mediated by the interaction of the system with the bath. At zero temperature, as the system undergoes further quantum relaxation, it eventually reaches its global ground state, which is also a pure state, the latter being the reason for the increase in the quantum purity at later times.

The overall shape of the different curves being the same, the main deviations between purities related to different coupling strengths appear in the magnitude of the corresponding minima, the specific points in time at which such minima are reached as well as in their asymptotic behavior. In the case of zero coupling, the system reaches its purity minimum around 5 ps after the dynamics have started, and its value at this particular instant is the largest among the minima appearing for all the other coupling strengths studied here. This indicates that the absence of inter-molecular vibrational coupling between the independent effective one-dimensional systems initially speeds up the loss of quantum coherence reaching much faster its minimum, but also the return to the pure ground state at later times. On the other hand, the purities associated with non-vanishing coupling strengths reach their minimum at much later times, around 10 ps. The presence of the inter-vibrational coupling between the effective onedimensional systems leads to a slower dynamics, as the two degrees of freedom internally exchange energy, delaying the energy release into the bath modes. This effect can be visualized both in the less pronounced slope of the purity when the inter-coupling constant is different from zero, as well as in the final value of the purity depicted in the graphic, as compared with the zero coupling case. Moreover, all of non-zero coupling purities have approximately the same value at their respective minimum, which corresponds to the most strongly mixed state. As the coupling strength increases, there are no major differences between the overall behavior of this magnitude along the entire simulation, except for the case of zero coupling. Regarding the comparison between the full and the factorizable methods, it can be observed from this figure that, as the other observables calculated in this work, both are very similar to each other for all cases. The small difference between the two approaches increases about linearly with coupling strength, apart from the case a = 0.1. Likewise, the asymmetric initial condition (2, 0) is responsible for the emergence of the structured pattern in the case of the fully coupled basis representation, resulting in small oscillations in the purity. In the small potential coupling regime, the incoherent relaxation due to the systembath interaction dominates, resulting in two different relaxation scenarios depending on the choice of the initial wave function representation. When the factorizable representation is selected, the initial wave-packet (2, 0) is an eigenstate of the effective Hamiltonian, Heff , and the initial density matrix is not only diagonal, but contains only one term different from zero, corresponding to the population of this single eigenstate. This feature leads to a nearly one-dimensional relaxation of the q1 degree of freedom at short times for systems where the characteristic environment-induced relaxation time is faster than or on the same order as IVR. Thus, the early system evolution is roughly dominated by the transitions (2, 0) → (1, 0) → (0, 0) essentially due to the system-bath coupling. In the factorizable ansatz, the resulting purity exhibits a smooth profile due to the vanishing effective inter-mode coupling. On the other hand, when the full spectral basis representation is chosen, this particular initial wave function is, by construction, a linear coherent superposition of

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

L. Uranga-Piña and J. C. Tremblay

dm(2,0) |m. Conversely, the

full

fact 2

√ < |p(t)n - p(t)n | >T

ij

imal spectral basis, two factorizable states √ combine to lift the degeneracy: |g = {(2, 0) + (0, 2)}/ 2 and |s = {(2, 0) √ − (0, 2)}/ 2, where the symmetric state |g is more stable. In general, since each spectral state exhibits both (2, 0) and (0, 2) character, the initial reduced density matrix associated with the factorizable initial condition is full and the states are coupled. These two-dimensional states implicitly carry the intermolecular vibrational coupling, thus leading to a somewhat distinct relaxation scenario where the two-dimensional wave function relaxes as a whole leading to natural oscillations in the purity. This reflects the coherent movement of the wave packet between the symmetrically equivalent factorizable states. These oscillations are damped after at longer times (t > 5 ps) due to the loss of “internal” coherence caused by the system-bath interaction. It is interesting to recognize that the frequency of this purity oscillation coincides with the corresponding oscillations of the populations. As the coupling strength increases (a = 0.3) an intermediate regime is reached, where the system-bath interaction and the inter-vibrational coupling are of the same order of magnitude and, although the purity undergoes the same behavior, its oscillation amplitude are almost vanishing from the very beginning. For larger inter-mode coupling constants, a, a regime where IVR is faster than the characteristic system-bath interaction time is reached. In this case the corresponding purities are almost flat (the oscillation frequencies are greater and the damping times of their amplitudes are smaller) and better coincides in the early stage with the purities obtained within the factorizable representation. Nevertheless, at longer times, the contribution of the loss of “internal” coherence (i.e., the IVR) within the full representation adds up, slightly slowing down the relaxation dynamics as compared to the pseudo-spectral representation. To gain more insight into the differences between the two theoretical approaches described above, the distances between two of the observables integrated over time as a function of the coupling strengths are shown in Fig. 8. The top panel depicts the average distances between the populations corresponding to each state along both modes, while the bottom panel depicts the distance between purities. Note that, for the case a = 0.1, the distances between the populations of the nq = 2 states for the two approaches does not exhibit i the same symmetry as for the other quantum numbers and the other coupling strengths. We attribute this difference to the asymmetric initial condition at small coupling strengths, as discussed above. In the coupled basis, the initial (2, 0) state is a linear superposition of eigenstates which includes large part of the IVR between states (2, 0) and (0, 2). This is also the case for the relaxation operators, which also include IVR. In the factorizable basis, only the potential contribution to the IVR is included and the relaxation operators remain decoupled. For weak potential-induced inter-mode coupling, it appears that system-bath induced coupling included in the dissipative operators can play a non-negligible role in

nq = 0 1 nq = 1 1 nq = 2 1 nq = 0 2 nq = 1 2 nq = 2

0.03 0.02

2

0.01 0 0.04

| >T

fully coupled states can be expressed as a linear combination of pseudo-spectral states, |m = cijm (i, j ). In a min-

0.04

fact 2

m

- γ(t)



full

two-dimensional states (2, 0) =

J. Chem. Phys. 141, 074703 (2014)

√ < |γ (t)

074703-12

0.03 0.02 0.01 0

0

0.1

0.2

0.3

0.4 0.5 0.6 coupling strength (a)

0.7

0.8

0.9

1

FIG. 8. Overall distances between observables computed using both approaches, for each inter-mode coupling strength a. (Top panel) Average distances between populations corresponding to each effective one-dimensional state along both modes. (Bottom panel) Distance between purities. The relative differences in all cases do not exceed the 4%.

the proper microscopic description of relaxation dynamics in the early stages. Once the system has lost a full vibrational quantum incoherently, this particular resonant condition is lost and the dynamics follow the same trend in both bases. As one increases the potential coupling, a, the energy transfer is faster and the resonant condition is reached earlier. This means that both the density matrices and the populations approach each other at early stages, thus causing the asymmetry feature to disappear more quickly. Also, the same mechanism explains why at intermediate times the distances between populations and purities corresponding to both methodologies are very similar. For larger coupling strengths, the energy transfer between equivalent states in the coupled basis representation is faster and accumulates, leading to a slightly slower relaxation dynamics at later times that cannot be fully described within the factorizable approach. Nonetheless, in both cases, the relative difference between the observables computed with both the full and the factorizable methods lies bellow 4%, thereby defining an upper limit for the validity of the factorizable ansatz, where it is possible to describe the dissipative quantum dynamics of a system interacting with its environment by explicitly constructing the corresponding multidimensional Lindblad operators as a tensor product of one-dimensional dissipative operators, and transferring the coupling between the modes to the remaining term of the commutator between the coupling potential and the density matrix in the Liouville-von Neumann equation (see Eq. (16)).

VI. CONCLUSIONS

We have analyzed the accuracy of approximating the dissipative operators appearing in the Liouville-von Neumann equation by a tensor product of effective one-dimensional operators carrying on most of the potential coupling between the different degrees of freedom of the system. The

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-13

L. Uranga-Piña and J. C. Tremblay

residual intra-molecular vibrational coupling term is taken into account by the commutator of the residual potential with the RDM in Eq. (16). The new approach was applied to the study of a two-dimensional system consisting in two CO molecules approaching a Cu(100) surface (Fig. 1). A comparison between the resulting dynamics using the factorizable method proposed in this work and the propagation of the reduced density matrix using fully coupled multidimensional Lindblad operators has been performed. The proposed method yields results in close agreement with the multidimensional propagation scheme, thereby confirming the validity of the factorization ansatz. The advantage of the form of the dissipative operators presented in this work over the standard spectral approach, lies in the possibility of treating larger systems that would be otherwise inaccessible. This advantage comes essentially from the feasibility of computing the multidimensional transition rates as a sum of one-dimensional rates, thus only requiring the diagonalization of N effective one-dimensional Hamiltonians instead of the much more computationally demanding diagonalization of the full N-dimensional Hamiltonian. In order to properly describe the dynamics of the system, the one dimensional dissipative operators contained inside the multidimensional Lindblad operator are designed to effectively account for most of the intra-molecular vibrational coupling. In the present study, such goal is achieved by a careful choice of the basis set used to represent these operators, that is, by employing a contracted grid basis generated via a PO-DVR transformation. The resulting partition of the total Hamiltonian of the system tends to flatten the attractive part of the coupling, so the residual interaction gets smaller around the active region where the dynamics takes place. It is worth to notice, that the resulting scheme is a factorizable, pseudo-spectral ansatz employing fully non-local relaxation operators, rather than a grid-based ansatz. The PO-DVR grid is used to define a good factorizable, pseudo-spectral representation that captures most of the intermode coupling. The one-dimensional populations computed within the two approaches described above, and starting from an initial wave-packet which is taken to be an eigenstate, (2, 0), of the effective one-dimensional Hamiltonians, agree remarkably well, even for the strongest inter-mode coupling. Regarding the case of zero coupling, they exactly coincide, owing to the fact that both approaches concur at this limit. The relaxation dynamics for zero inter-vibrational coupling occurs only along the q1 axis, while decoherence is faster than for non-vanishing inter-mode couplings. For a > 0.0, the relaxation mechanism in the factorizable picture consists initially in population going back and forth from state (2, 0) to state (0, 2) via state (1, 1) during the first 3 ps, while slowly relaxing to the asymptotic (0, 0) state. The oscillatory pattern exhibited by these one-dimensional populations can be regarded as Rabi oscillations, similar to the problem of a three level system subject to an external oscillatory field. Within this picture, the role of the external oscillatory field is played by the average external oscillatory force exerted by the remaining degrees of freedom. The greater the inter-mode coupling, the larger the Rabi frequency observed in the population patterns.

J. Chem. Phys. 141, 074703 (2014)

The time-dependent probability density distributions for the whole range of inter-vibrational couplings considered in this work are remarkably similar for both the factorization ansatz and the full spectral representation. The main features of the relaxation dynamics of the system, e.g., the inversion of population with respect to the initial distribution at some points in time, are well reproduced. Also, the overall dynamics exhibit periodic breathings of the wave packet along each of the coordinates for the first 11 ps, eventually relaxing to the global ground state (0, 0). Regarding the quantum purity we have observed that there are minor differences between both descriptions for the whole range of potential couplings studied. Since the initial wave-packet coincides with one of the eigenstates of the factorizable Hamiltonian, the purity for small inter-mode couplings initially exhibits some small oscillatory features in the fully coupled spectral basis that are absent in the factorizable ansatz. We attribute this behavior to the fast relaxation along a single mode as compared to the inter-mode coupling time within the Lindblad operators. Owing to the zero temperature regime, the system purity is fully recovered in the asymptotic limit in all cases, which occurs faster in absence of IVR. The presence of the vibrational coupling between the modes slows down the process of decoherence and the recovering of the asymptotic coherent character, since the energy exchange between the internal modes delays the energy release into the bath modes. Based on the results summarized above, it can be concluded that the proposed approach of representing the dissipative Lindblad operators as a tensor product of onedimensional operators, combined with a suitable potentialoptimized basis representation, provides a proper description of the population dynamics and probability distributions evolution in quantum dissipative systems. The two-dimensional system addressed here incorporates the essential features of the dynamics of adsorbates at surfaces. Hence, the close correspondence between the results of the present model and those of standard techniques for the integration of the Liouville-von Neumann equation, is a strong indication that most of the intermode coupling upon relaxation is mediated by IVR and not by the coupling between relaxation operators, even at large coupling strengths. Since the former can be included exactly in a factorizable ansatz, the conclusions presented here are likely to hold for dissipative dynamics in similar systems. Compared to the use of coupled multidimensional spectral Lindblad operators, the present scheme presents the advantage of enabling a marked reduction of the grid size and a favourable scaling of the computational cost with respect to the increase of the dimensionality of the system. In the present study, the vibrational degrees of freedom other than the CO-surface stretch mode, have been disregarded. As an outlook, we plan to consider other degrees of freedom of the adsorbate-surface system, such as frustrated translations and rotations. The inclusion of the remaining vibrational modes is a necessary step for an accurate description of the true system dynamics, and it constitutes the natural extension of the present algorithm to treat higher dimensional systems.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-14

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

ACKNOWLEDGMENTS

full

fact. (t = 0.96 ps) q 2 [a.u.]

5.2

See the snapshots of the density evolution in Fig. 9.

5.0 4.8 4.6

(t = 3.36 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

(t = 4.80 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

5.2

(t = 10.80 ps) q 2 [a.u.]

(t = 10.80 ps) q 2 [a.u.]

(t = 5.52 ps) q 2 [a.u.]

(t = 2.88 ps) q 2 [a.u.]

(t = 0.96 ps) q 2 [a.u.]

The authors gratefully acknowledge the funding of the Deutsche Forschungsgemeinschaft (DFG) through the Emmy-Noether program, Project No. Tr1109/2-1.

APPENDIX: EVOLUTION OF THE DENSITY FOR VARIOUS COUPLING STRENGTHS

5.0 4.8 4.6 4.6 4.8 5.0 5.2 q 1 [a.u.]

5.2

full

fact.

4.6 4.8 5.0 5.2 q 1 [a.u.]

4.6 4.8 5.0 5.2 q 1 [a.u.]

5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

4.6 4.8 5.0 5.2 q 1 [a.u.]

full

fact. (t = 0.24 ps) q 2 [a.u.]

5.2

(b)

5.0 4.8 4.6

(t = 2.16 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

(t = 4.08 ps) q 2 [a.u.]

5.2 5.0 4.8 4.6

5.2

(t = 10.80 ps) q 2 [a.u.]

(t = 10.80 ps) q 2 [a.u.]

(t = 5.76 ps) q 2 [a.u.]

(t = 2.64 ps) q 2 [a.u.]

(t = 0.24 ps) q 2 [a.u.]

(a)

5.0 4.8 4.6 4.6 4.8 5.0 5.2 q 1 [a.u.]

4.6 4.8 5.0 5.2 q 1 [a.u.]

(c)

5.2

full

fact.

4.6 4.8 5.0 5.2 q 1 [a.u.]

4.6 4.8 5.0 5.2 q 1 [a.u.]

5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

5.2 5.0 4.8 4.6

(d)

FIG. 9. Snapshots of the density evolution for an initial condition (2, 0). Weak coupling (a) (a = 0.1) and (b) (a = 0.3). Strong coupling (c) (a = 0.8) and (d) (a = 1.0)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

074703-15

L. Uranga-Piña and J. C. Tremblay

J. Chem. Phys. 141, 074703 (2014)

1 B.

30 N.

2 J.

31 J.

Hammer, Adv. Catal. 45, 71 (2000). Greeley, J. K. Norskov, and M. Mavrikakis, Annu. Rev. Phys. Chem. 53, 319 (2002). 3 K. Mayrhofer, V. Juhart, K. Hartl, M. Hanzlik, and M. Arenz, Angew. Chem., Int. Ed. 48, 3529 (2009). 4 U. Panzer and H. P. Schreiber, Macromolecules 25, 3633 (1992). 5 J. B. Donnet, S. J. Park, and H. Balard, Chromatographia 31, 434 (1991). 6 S. J. Lee and M. Moskovits, Nano Lett. 11, 145 (2011). 7 J. C. Tully, Annu. Rev. Phys. Chem. 51, 153 (2000). 8 N. Camillone III, K. A. Khan, P. J. Lasky, L. Wu, J. E. Moryl, and R. M. Osgood Jr., J. Chem. Phys. 109, 8045 (1998). 9 M. Bonn, S. Funk, C. Hess, D. N. Denzler, C. Stampfl, M. Scheffler, M. Wolf, and G. Ertl, Science 285, 1042 (1999). 10 N. Koumura, R. W. J. Zijlstra, R. A. van Delden, N. Harada, and B. L. Feringa, Nature 401, 152 (1999). 11 C. J. Hirschmugl and G. P. Williams, Phys. Rev. Lett. 65, 480 (1990). 12 F. Hofmann and J. P. Toennies, Chem. Rev. 96, 1307 (1996). 13 A. P. Graham, F. Hofmann, J. P. Toennies, G. P. Williams, C. J. Hirschmugl, and J. Ellis, J. Chem. Phys. 108, 7825 (1998). 14 M. Bonn, C. Hess, S. Funk, J. H. Miners, B. N. J. Persson, M. Wolf, and G. Ertl, Phys. Rev. Lett. 84, 4653 (2000). 15 E. H. G. Backus, A. Eichler, A. W. Kleyn, and M. Bonn, Science 310, 1790 (2005). 16 K. D. Dobbs and D. J. Doren, J. Chem. Phys. 99, 10041 (1993). 17 P. Hu, D. A. King, S. Crampin, M. H. Lee, and M. C. Payne, J. Chem. Phys. 107, 8103 (1997). 18 S. Carter, S. J. Culik, and J. M. Bowman, J. Chem. Phys. 107, 10458 (1997). 19 H. Guo and R. Chen, J. Chem. Phys. 110, 6626 (1999). 20 M. Berman, R. Kosloff, and H. Tal-Ezer, Phys. A: Math. Gen. 25, 1283 (1992). 21 I. Andrianov and P. Saalfrank, Chem. Phys. Lett. 350, 191 (2001). 22 I. Andrianov and P. Saalfrank, J. Chem. Phys. 124, 034710 (2006). 23 B. Gergen, H. Nienhaus, W. H. Weinberg, and E. W. McFarland, Science 294, 2521 (2001). 24 J. Y. Park and G. A. Somorjai, Chem. Phys. Chem. 7, 1409 (2006). 25 S. Gao, M. Persson, and B. I. Lundqvist, Solid State Commun. 84, 271 (1992). 26 K. Golibrzuch, A. Kandratsenka, I. Rahinov, R. Cooper, D. J. Auerbach, A. M. Wodtke, and C. Bartels, J. Phys. Chem. A 117, 7091 (2013). 27 B. Mildner, E. Hasselbrink, and D. Diesing, Chem. Phys. Lett. 432, 133 (2006). 28 H. Kasai and A. Okiji, Surf. Sci. Lett. 225, L33 (1990). 29 S. Gao, M. Persson, and B. Lundqvist, J. Electron. Spectrosc. Relat. Phenom. 64-65, 665 (1993).

Shenvi, S. Roy, and J. C. Tully, J. Chem. Phys. 130, 174107 (2009). C. Tully, M. A. Gomez, and M. Head-Gordon, J. Vac. Sci. Technol., A 11, 1914 (1993). 32 M. Nest and P. Saalfrank, J. Chem. Phys. 116, 7189 (2002). 33 P. Saalfrank, Chem. Rev. 106, 4116 (2006). 34 G. Lindblad, Commun. Math. Phys. 48, 119 (1976). 35 K. Blum, Density Matrix Theory and Applications (Plenum Press, New York, 1996). 36 H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, New York, 2002). 37 M. Razavy, Classical and Quantum Dissipative Systems (Imperial College Press, London, 2005). 38 A. S. Leathers and D. A. Micha, J. Phys. Chem. A 110, 749 (2006). 39 A. S. Leathers, D. A. Micha, and D. S. Kilin, J. Chem. Phys. 131, 144106 (2009). 40 C. Meier and D. Tannor, J. Chem. Phys. 111, 3365 (1999). 41 U. Kleinekathöfer, J. Chem. Phys. 121, 2505 (2004). 42 M. Head-Gordon and J. C. Tully, J. Chem. Phys. 96, 3939 (1992). 43 J. T. Kindt, J. C. Tully, M. Head-Gordon, and M. A. Gomez, J. Chem. Phys. 109, 3629 (1998). 44 B. Hellsing and M. Persson, Phys. Scr. 29, 360 (1984). 45 J. C. Tremblay and P. Saalfrank, J. Chem. Phys. 131, 084716 (2009). 46 J. C. Tremblay, S. Monturet, and P. Saalfrank, Phys. Rev. B 81, 125408 (2010). 47 J. C. Tremblay, J. Chem. Phys. 138, 244106 (2013). 48 M. H. Beck, A. Jackle, G. Worth, and H. D. Meyer, Phys. Rep. 324, 1 (2000). 49 D. V. Shalashilin and M. S. Child, J. Chem. Phys. 115, 5367 (2001). 50 P. Saalfrank, G. Boendgen, K. Finger, and L. Pesce, Chem. Phys. 251, 51 (2000). 51 K. Finger and P. Saalfrank, Chem. Phys. Lett. 268, 291 (1997). 52 S. Li and H. Guo, J. Chem. Phys. 117, 4499 (2002). 53 V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976). 54 A. Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer and Reactions in Condensed Molecular Systems (Oxford University Press, New York, 2006). 55 A. Rajagopal, Phys. Lett. A 246, 237 (1998). 56 R. Marquardt, F. Cuvelier, R. A. Olsen, E. J. Baerends, J. C. Tremblay, and P. Saalfrank, J. Chem. Phys. 132, 074108 (2010). 57 J. C. Tremblay, G. Füchsel, and P. Saalfrank, Phys. Rev. B 86, 045438 (2012). 58 J. Echave and D. C. Clary, Chem. Phys. Lett. 190, 225 (1992). 59 B. W. Shore, Manipulating Quantum Structures Using Laser Pulses (Cambridge University Press, New York, 2011). 60 T. A. Laine and S. Stenholm, Phys. Rev. A 53, 2501 (1996).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Thu, 27 Nov 2014 17:56:35

Relaxation dynamics in quantum dissipative systems: the microscopic effect of intramolecular vibrational energy redistribution.

We investigate the effect of inter-mode coupling on the vibrational relaxation dynamics of molecules in weak dissipative environments. The simulations...
2MB Sizes 0 Downloads 3 Views