Quantum Stochastic Trajectories: The Fokker–Planck–Bohm Equation

Feb 21, 2018 - To overcome such a drawback, we derive a stochastic representation of the quantum molecular trajectory, through projection operator tec...
0 downloads 7 Views 1MB Size
Subscriber access provided by UNIV OF DURHAM

Article

Quantum Stochastic Trajectories: the Fokker-PlanckBohm Equation Driven by the Reduced Density Matrix Francesco Avanzini, and Giorgio J. Moro J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b11943 • Publication Date (Web): 21 Feb 2018 Downloaded from http://pubs.acs.org on March 1, 2018

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

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

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

The Journal of Physical Chemistry

Quantum Stochastic Trajectories: the Fokker-Planck-Bohm Equation Driven by the Reduced Density Matrix Francesco Avanzini∗ and Giorgio J. Moro∗ Dipartimento di Scienze Chimiche, Università di Padova, via Marzolo 1, 35131 Padova, Italy E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment

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

Abstract The quantum molecular trajectory is the deterministic trajectory, arising from the Bohm theory, that describes the instantaneous positions of the nuclei of molecules by assuring the agreement with the predictions of Quantum Mechanics. Therefore, it provides the suitable framework for representing the geometry and the motions of molecules without neglecting their quantum nature. However, the quantum molecular trajectory is extremely demanding from the computational point of view, and this strongly limits its applications. To overcome such a drawback, we derive a stochastic representation of the quantum molecular trajectory, through projection operator techniques, for the degrees of freedom of an open quantum system. The resulting Fokker-Planck operator is parametrically dependent upon the reduced density matrix of the open system. Because of the pilot role played by the reduced density matrix, this stochastic approach is able to represent accurately the main features of the open system motions both at equilibrium and out of equilibrium with the environment. To verify this procedure, the predictions of the stochastic and deterministic representation are compared for a model system of six interacting harmonic oscillators, where one oscillator is taken as the open quantum system of interest. The undeniable advantage of the stochastic approach is that of providing a simplified and self-contained representation of the dynamics of the open system coordinates. Furthermore, it can be employed to study the out of equilibrium dynamics and the relaxation of quantum molecular motions during photoinduced processes, like photoinduced conformational changes and proton transfers.

1

Introduction

The same idea of molecules, as objects characterized by a time dependent geometry describing their motions in the three dimensional space, implies the assignment to their nuclei of well defined spatial positions evolving in time. This is coherent with the use of Classical Mechanics to describe Molecular Dynamics, since it allows the univocal identification of instantaneous nuclear positions. In this framework the conformational dynamics of biomolecules is 2

ACS Paragon Plus Environment

Page 2 of 42

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

The Journal of Physical Chemistry

often analyzed in order to recognize their biological role. Also chemical reactions are usually interpreted as the movements of particular molecular fragments, like the proton transfer in simple acid-base reactions. On the other hand, the nuclei are quantum particles like the electrons. Especially for light nuclei, nuclear quantum effects can not be neglected and different methods have been developed to take them into account, like the ring-polymer molecular dynamics, 1,2 the semiclassical initial value representation 3 and the centroid molecular dynamics. 4,5 For example, the pH of water would be 8.5 instead of 7, if the protons are treated classically. 6 However, the classical and the quantum representations are not mutually consistent: their predictions do not correspond and, in addition, Quantum Mechanics does not allow the univocal identification of the spatial coordinates for any kind of particle, even the nuclei. Recently we have proposed what has been called the quantum molecular trajectory to determine the time dependent particle positions without contradicting at the same time the predictions of Quantum Mechanics. 7 The later are recovered in terms of statistical properties of the quantum molecular trajectory. For example, the expectation values are formally related to the time averages of the corresponding physical observables along the trajectory. This leads to a rigorous representation of the time dependent molecular geometry within a fully quantum framework. In Ref. 7 we have formally proven that the quantum molecular trajectory can be identified with a single Bohm trajectory. The Bohm theory is a formulation of Quantum Mechanics that characterizes the state of a quantum system according to both the wave function, as in the standard formulation, and the time dependent configuration. 8–10 The configuration defines the spatial position of all the particles like in Classical Mechanics and it evolves in time by drawing a quantum continuous trajectory, i.e., the Bohm trajectory. In the original formulation by Bohm, 8,9 the equivalence with Quantum Mechanics is recovered by employing the ensemble of all possible configurations evolving along a swarm of trajectories that collectively reproduce the time dependence of the wave function. This equivalence has been extensively exploited for developing innovative algorithms to calculate the

3

ACS Paragon Plus Environment

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

wave function dynamics 11,12 and to derive semi-classical approximations. 13,14 As A. Figalli and coworkers stated, 15 despite Bohm theory “as particle trajectories remains controversial from the physics point of view, its mathematical foundation is solid”: the Bohm trajectory is well defined for all times. 16,17 In our approach we employ precisely a single Bohm trajectory, in opposition to the statistical ensemble, for identifying the nuclear positions. The undeniable benefit of the quantum molecular trajectory is that of dealing with a trajectory for representing the molecular motions that does not require any correction in order to preserve the coherence with Quantum Mechanics. However, this theoretical advantage is strongly limited by the computational cost of solving the dynamical equations. On the one hand, the wave function of the isolated quantum system must be known at any time through the solution of the Schrödinger equation. On the other hand, given the wave function determining the velocity of the particles according to the Bohm equation, the quantum trajectory has to be determined by solving ordinary differential equations that are characterized by a rather stiff behaviour. As a matter of fact, the quantum molecular trajectory can be calculated only for isolated systems composed of few degrees of freedom. In this work, we intend to overcome these computational drawbacks, by deriving a stochastic representation of the quantum molecular trajectory. Stochastic methods are shared by several fields of Physical Chemistry. 18 They can be employed for example in the study of the kinetics of chemical reactions by considering explicitly the fluctuations on the number of molecules of the components. 19–22 Recently, an analysis of two-states kinetics that takes into account the noise effects on the measured concentrations has been proposed. 23 Furthermore, it is widely known the importance of stochastic methodologies for the investigation of the dynamics of macromolecules of biological interest, such as proteins. 24,25 The correlation between overall motion and the internal dynamics of proteins has been discussed recently. 26–28 Our aim concerns the stochastic representation of the Bohm coordinate dynamics and, therefore, from the methodological point of view, it is completely different from the dissipative or stochastic extension of the Schrödinger equation. 29–33

4

ACS Paragon Plus Environment

Page 4 of 42

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

The Journal of Physical Chemistry

In particular, we focus on the degrees of freedom of an open quantum system of interest interacting with a much larger environment instead of all the set of degrees of freedom of a closed (isolated) system. For example, when a conformational change occurs, the degrees of freedom characterizing the molecular geometry are more relevant in describing the process that the others. By employing projection operator techniques, like Mori-Zwanzig procedures in Classical Mechanics, 34–37 we derive the stochastic equations for the set of relevant degrees of freedom from the Liouville representation of the Schrödinger-Bohm dynamics. 7 By considering the open system coordinates as the only relevant variables, in Ref. 38 we have derived a Smoluchowski type equation, the so-called Smoluchowski-Bohm equation, for the distribution of the Bohm coordinates of the open system. In this framework the dynamics of the wave function does not appear explicitly besides its role as a source of the noise. We have verified that Smoluchowski-Bohm equation supplies a satisfactory approximation of the quantum molecular trajectory picture when the open system is in equilibrium with the environment. On the other hand, the wave function evolution plays a fundamental role when the open system is out of equilibrium and in this case one needs a more general representation. In this work, we derive a new stochastic representation that includes explicitly also the wave function dynamics. To this aim we use a different projection operator that leads to a Fokker-Planck operator acting on the density distribution of the Bohm coordinates of the open system, and which is parametrically dependent on the corresponding reduced density matrix. We refer to the resulting stochastic representation as the Fokker-Planck-Bohm equation driven by reduced density matrix, whose ingredients are supplied by standard Quantum Mechanics. For example, the equilibrium distribution of the stochastic process is replaced by the wave function’s squared modulus integrated on the degrees of freedom of the environment. In order to validate our stochastic representation, we compare the predictions of the stochastic equations with those of the deterministic quantum molecular trajectory of the isolated system. To this purpose, we employ a model system of six interacting harmonic oscillators initially out of equilibrium, where one oscillator plays the role of the open system

5

ACS Paragon Plus Environment

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

of interest while the others represent the environment. The benefits of our stochastic procedure based on the Fokker-Planck-Bohm equation concern both methodological and computational issues. From the methodological point of view, it represents a convenient quantum framework for describing the molecular motions in terms of nuclear quantum trajectories like the deterministic quantum molecular trajectory for the isolated system. 7 From the computational point of view, the Fokker-Planck-Bohm equation supplies a reduced representation of the dynamics by considering only the Bohm coordinates of the open system with the reduced density matrix playing the role of the pilot filed, very much like the wave function in the Bohm equation for the isolated system. In this way our method preserves the low computational cost that is typical of the stochastic methods. These features should make our description based on the Fokker-Planck-Bohm equation very useful for characterizing the molecular dynamics. Furthermore, the capability of tackling the out of equilibrium dynamics should make it suitable to investigate complex dynamical processes like photoinduced conformational changes (as the photoisomerization of azobenzene) or photoinduced proton transfers. The manuscript is organised as follow. In Section 2 we summarize the theoretical background. In particular, the Schrödinger-Bohm dynamics is considered in subsection 2.1, while the statistical characterization of the quantum molecular trajectory is introduced in subsection 2.2. In the subsections 2.3, 2.4, 2.5 the stochastic procedure, the Fokker-Planck-Bohm equation and the role of the reduced density matrix are examined, respectively. We present the model system of six interacting harmonic oscillators in Section 3, while the comparison between the stochastic and the deterministic description is made in Section 4. Conclusion are drawn in Section 5.

6

ACS Paragon Plus Environment

Page 6 of 42

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

The Journal of Physical Chemistry

2

Theory

The objective of our theoretical analysis is the stochastic description of the Bohm coordinates for the open quantum system interacting with the environment. To this aim, the formal tools already introduced in Ref. 7 and in Ref. 38 are employed in order to derive the Fokker-Planck equation for the coordinates of the open quantum system. In the following three subsections, these tools are recalled before to examine the final Fokker-Planck description.

2.1

Schrödinger-Bohm dynamics

Bohm dynamics 8,9 can be understood as a formal extension of the conventional Schröedinger dynamics, 15 which supplies the time dependent configuration Q(t) of a closed (isolated) quantum system on the basis of its quantum pure state |Ψ(t)i having a pilot role. For the quantum system with n cartesian degrees of freedom, the notation Q(t) = (Q1 (t), Q2 (t), · · · , Qn (t)) will be used for the instantaneous configuration, where Qk (t) is the coordinate of the k-th degrees of freedom (k = 1, 2, · · · , n). A deterministic description of the overall dynamics is provided by the following set of differential equations ) ( ∗  ∂ Ψ(q, t) Ψ (q, t)  ~ d ∂qk   Q (t) = Im k  dt 2 mk |Ψ(q, t)| q=Q(t)    ∂ ı  |Ψ(t)i = − H|Ψ(t)i ˆ ∂t ~

(1)

ˆ is the Hamiltonian where mk is the mass corresponding to the k-th degree of freedom, H operator and the wave function Ψ(q, t) supplies the coordinate representation of the pure state, Ψ(q, t) = hq|Ψ(t)i, where the coordinates q are specified as q = (q1 , q2 , · · · , qn ). It is evident that the evolution of Bohm coordinates is driven by Ψ(q, t) without any back effect. The correspondence between the statistical properties of Bohm coordinates and the quantum pure state |Ψ(t)i is recovered if one considers a swarm of trajectories Q(t) initially distributed as |Ψ(q, 0)|2 . As shown by Bohm in his original contributions, 8,9 the distribution generated by such a swarm of trajectories on the configuration space at an arbitrary time t is equivalent to the coordinate distribution |Ψ(q, t)|2 determined by the wave function. This 7

ACS Paragon Plus Environment

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

Page 8 of 42

equivalence is exploited in numerical procedures for the calculation of the time dependent wave function Ψ(q, t) from a suitable ensemble of trajectories. 11,39 Our methodological point of view, however, is rather different: the Bohm coordinates are employed to characterize the quantum molecular trajectory as long as they allow the assignment of well defined positions to the particles of the system, the nuclei of a molecule for instance. Correspondingly the single Bohm trajectory Q(t), in principle to be evaluated for a given time dependent pure state |Ψ(t)i and initial configuration Q(0), represents the main objective of our analysis. 40 As shown in Ref. 7, the correspondence with the quantum description provided by the wave function is recovered from the distribution on the configuration space generated by a single Bohm trajectory Q(t) during its evolution.

2.2

Liouville representation of the deterministic dynamics

Let us parameterize the time dependent pure state |Ψ(t)i according to the eigenstates {|φj i} ˆ j i = Ej |φj i for j = 1, 2, . . . . We and to the eigenvalues {Ej } of the Hamiltonian: H|φ assume that the energy eigenvalues are rationally independent 41 and that the wave function has components along the first N eigenstates only. We do not impose any a priori constraint on the number of eigenstates and, therefore, the following analysis can be applied to any quantum state characterized by a finite number N of its eigenstate components. Thus the pure state can be specified as |Ψ(t)i =

N X p Pj e−ıAj (t) |φj i,

(2)

j=1

with constant and normalized populations P = (P1 , P2 , · · · , PN ) of the eigenstates,

PN

j=1

Pj =

1, and with linearly time dependent phases A(t) = (A1 (t), A2 (t), · · · , AN (t)) where Aj (t) = Aj (0) + Ej t/~. Correspondingly, for a given set P of populations, the dynamical state X(t) evolving deterministically according to eq. (1) can be specified at each instant t by means of the set of the n Bohm coordinates Q(t) together with the set of the N quantum phases

8

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

A(t), X(t) = (Q(t), A(t)) .

(3)

We introduce also the invariant manifold ΩP for the ensemble of possible dynamical states X(t), with reference to a specific set P of populations. 7 A generic point of ΩP will be denoted as x = (q, α). The invariant manifold ΩP represents the natural support for the distribution on the dynamical states, more specifically the time dependent probability density ρΩP (x, t) to be normalized by integration over ΩP . As a consequence of the deterministic dynamics of X(t), such a distribution follows the Liouville equation 7,42 ∂ ΩP ρ (x, t) = −∇P · ΛP (x)ρΩP (x, t), ∂t

(4)

where ∇P is the gradient operator for the x variables, while the n + N components of the vector field ΛP (x) (or velocity field) are given as ~ ΛPk (q, α) = Im mk

(

ψP∗ (q, α) ∂q∂k ψP (q, α) |ψP (q, α)|2

)

,

ΛPn+j = Ej /~,

(5) (6)

with k = 1, 2, . . . , n and j = 1, 2, . . . , N . In the previous relation the symbol ψP (q, α) is used to denote the wave function Ψ(q, t) = hq|Ψ(t)i recovered from eq. (2) by means of the coordinate representation of the Hamiltonian eigenstates φj (q) = hq|φj i, with an implicit time dependence brought by the phases:

ψP (x) x=(q,α) = ψP (q, α) :=

N X p Pj e−ıαj φj (q),

(7)

j=1

so that the pure state wave function at a given time t is recovered by inserting the corresponding values of the phases: Ψ(q, t) = ψP (q, α)|α=A(t) . In eq. (4), as well as in the next equations, the dot symbol is employed to denote the scalar product between the components of the generic gradient operator ∇ and of the generic velocity field Λ(•). Finally, under the condition that the system Hamiltonian has the ordinary mechanical structure including the standard kinetic energy operator and the coordinate dependent 9

ACS Paragon Plus Environment

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

Page 10 of 42

potential energy, and by taking into account the specific form of the velocity field ΛP (x), see eq. (5) and eq. (6), one derives the following equilibrium distribution as the stationary solution of the Liouville Equation (4), 7 P ρΩ eq (x) =

1 |ψP (x)|2 (2π)N

(8)

P with the 1/(2π)N factor ensuring the normalization of the equilibrium distribution ρΩ eq (x)

by integration on the coordinates q and the phases α.

2.3

Stochastic representation of the dynamics

The full solution of deterministic equations (1) appear unfeasible in realistic models of material systems because of the huge number of degrees of freedom. Thus, an approximate representation of the dynamics has necessarily to be considered by focusing on the observables of interest. In all generality we introduce a partition of the variables of eq. (3) specifying the dynamical state X(t) = (XR (t), XI (t))

(9)

according to the set XR (t) = {Xi (t)}i∈R of relevant variables and the complementary set XI (t) = {Xi (t)}i∈I of irrelevant variables which are not required to evaluate the properties of interest. By integrating the Liouville distribution ρΩP (x, t) on the irrelevant variables, R

ρ (xR , t) :=

Z

dxI ρΩP (x, t)|x=(xR ,xI ) ,

(10)

one gets directly the marginal distribution ρR (xR , t) on the relevant variables. The same P integration can be applied to the Liouville equilibrium distribution ρΩ eq (x) to obtain the re-

duced equilibrium distribution ρR eq (xR ). Furthermore, projection operator techniques allows the formal derivation of an independent equation for the evolution of the reduced distribution ρR (xR , t) which is physically justified in the presence of a suitable time scale separation leading to a Markovian dynamics for the relevant variables. 34–37 In Ref. 38 we have described the projection procedure applied to the Liouville equation for the Schrödinger-Bohm dynamics and here we recall only the final results. 10

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

Without introducing any particular choice for the relevant variables, the following FokkerPlanck equation is derived in all generality for the reduced distribution h ∂ R ρ (xR , t) = − ∇R · ΛR (xR )+ ∂t −1 i R R − ∇R · D ρ R ρR ρ (xR , t), eq (xR )∇ eq (xR )

(11)

where the diffusion matrix D acts on the gradient operator ∇R with respect to the relevant variables. The whole term between square brackets [. . . ] in eq. (11) is the Fokker-Planck operator that acts on the density distribution ρR (xR , t) on the relevant variables. Like with the stochastic treatments of classical dynamics, 18 the Fokker-Planck operator is separated into a streaming-like contribution, the first term at the r.h.s. of the previous equation, and a collisional-like contribution, the second term at the r.h.s. of the same equation. The former is driven by the reduced velocity field ΛR (xR ) which is simply the set of the component of the deterministic velocity field ΛP (x) for the relevant variable xR after the average on the irrelevant ones: 38 broadly speaking, ΛR (xR ) represents the average drift for the relevant variables xR . The later term takes into account the fluctuating contributions to the dynamics of the relevant variable and, consequently, it describes the relaxation processes, such as the decay of a generic distribution ρR (xR , t) towards the equilibrium distribution ρR eq (xR ). The collisional term is controlled by the diffusion matrix D, accounting for the time scale and the amplitude of the fluctuations, whose formal definition as supplied by the projection procedure, however, cannot be employed to estimate its values. Therefore the diffusion matrix has to be considered as a set of free parameters to be optimized on the basis of the predictions of the stochastic representation. It should be mentioned that in all generality the xR dependence of the diffusion matrix has to be considered but, as usually done with classical stochastic methods, we suppose that such a dependence is negligible. In the following, we will specialize the generic Fokker-Planck operator of eq. (11) in order to address the specific objectives of this work, but the physical meaning of the streaming-like contribution and of the collisional-like contribution is the same for all the examined forms of this type of Fokker-Planck equations. 11

ACS Paragon Plus Environment

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

Page 12 of 42

The theory of stochastic equations establishes an equivalence between the description supplied by the probability distribution ruled by the Fokker-Planck equation and the description based on fluctuating trajectories described by stochastic differential equations of Langevin type. 18 Given the Fokker-Planck equation (11), the corresponding Langevin equation for the stochastic evolution of the relevant variables is specified as 38 d XR (t) = ΛR (XR (t)) + dt −1 1/2 + ρR D ∇R ρ R ξ(t), eq (XR (t)) eq (XR (t)) + (2D)

(12)

where ξ(t) is a multidimensional white noise.

2.4

The Fokker-Planck-Bohm equation

In this subsection we specialize the previous stochastic formalism to the open quantum system in order to derive a convenient representation of the dynamics of its Bohm coordinates. Let us consider the closed quantum system, which follows the deterministic Schröedinger-Bohm dynamics, as composed by the open quantum system S of interest and the environment-bath B. Correspondingly the full set Q of Bohm coordinates is decomposed into the set QS = {Qk }k∈S of coordinates for the open quantum system S, and the coordinates QB = {Qk }k∈B for the environment B: Q = (QS , QB ). The objective of our stochastic treatment is a dynamical representation of the QS coordinates independently of the environment. Also the quantum pure state |Ψ(t)i of the closed system can be specialized for the open quantum system by means of the reduced density matrix σ ˆ (t) defined as the trace on the environment states of the pure state density matrix |Ψ(t)ihΨ(t)|, that is σ ˆ (t) := TrB {|Ψ(t)ihΨ(t)|}. Of course σ ˆ (t) is an operator acting on Hilbert space HS for the open quantum system. If the α-parametrized representation of the pure state is adopted according to eq. (7), the following phase dependent reduced density matrix is recovered σ ˆ (α) =

N X p

j,j 0 =1

Pj Pj 0 e−ı(αj −αj 0 ) TrB {|φj ihφj 0 |},

12

ACS Paragon Plus Environment

(13)

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

The Journal of Physical Chemistry

from which the reduced density matrix at a given time is evaluated by inserting the actual values of the phases, σ ˆ (t) = σ ˆ (α)|α=A(t) . Let us consider the particular situation of an open quantum system in equilibrium with the environment. In this case the reduced density matrix σ ˆ (t) displays small amplitude P ˆ = N fluctuations around its time average σ j=1 Pj TrB {|φj ihφj |} which well reproduces the

expectation values for the open quantum system. 43,44 Then, it is sufficient to identify the relevant variables of the projection procedure with the Bohm cooordinates of the open quantum system, XR (t) = QS (t) = {Qk (t)}k∈S , in order to characterize the stochastic trajectory of

these coordinates. Correspondingly, as shown in detail in Ref. 38, a purely diffusive motion described by the so-called Smoluchowski-Bohm equation is recovered. In the present work we consider the more general situation of a non-equilibrium state due the evolution of the the quantum pure state |Ψ(t)i from a particular initial condition |Ψ(0)i. For instance one might imagine, like in non-linear optical experiments, to apply a strong external field to the open quantum system in order to generate an initial non-equilibrium state. This leads to a reduced density matrix σ ˆ (t) with an explicit time dependence before ˆ . In this case the coordinates of the open reaching the equilibrium condition described by σ quantum system should experience significant effects of the pure state evolution and the projection procedure with the restriction of the relevant variables XR to QS is unable to capture these effects. One should take care of the dynamics of the quantum pure state in the projection procedure, and this is simply done by including also the phases A(t) amongst the relevant variables, XR (t) = (QS (t), A(t)) .

(14)

Thus, the irrelevant variables XI are identified with the Bohm coordinates QB of the environment. Of course effects of environment coordinates QB on the open system coordinates QS cannot be excluded since, within the deterministic Schröedinger-Bohm representation eq. (1), the velocity field acting on QS depends through the wave function on all the Bohm coordinates Q(t) = (QS (t), QB (t)). However, because of the strongly chaotic behaviour 13

ACS Paragon Plus Environment

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

Page 14 of 42

of Bohm coordinates, 45–48 one expects that the coupling with the environment coordinates should produce only correlations at very short times which are compatible with a Markovian representation for the variables of eq. (14). Given the choice of eq. (14) for the relevant variables, according to eq. (10) the reduced distribution is obtained by integration on the environment coordinates ρR (qS , α, t) :=

Z

dqB ρΩP (q, α, t)|q=(qS ,qB ) .

(15)

An analogous relation defines the reduced equilibrium distribution ρR eq (qS , α), which can be further specified in terms of the elements of reduced density matrix σ ˆ (α) if a generic orthonormal basis {|ηlS i} is available for HS ρR eq (qS , α) =

:=

Z

P dqB ρΩ eq (q, α)|q=(qS ,qB )

1 X σlS0 ,lS (α)ηl∗S (qS )ηlS0 (qS ), (2π)N 0

(16)

lS ,lS

where σlS ,lS0 (α) = hηlS0 |ˆ σ (α)|ηlS i and ηlS (qS ) is coordinate representation of |ηlS i, ηlS (qS ) = hqS |ηlS i. Furthermore the generic Fokker-Planck Equation (11) can be specialized for the relevant variables of eq. (14) to obtain X N Ej ∂ ∂ R ρ (qS , α, t) = − + ∇S · ΛS (qS , α)+ ∂t ~ ∂αj j=1  −1 R R S R S ρ (qS , α, t). − ∇ · D ρeq (qS , α)∇ ρeq (qS , α)

(17)

In this equation ∇S is the gradient operator with respect to Bohm coordinates qS of the open quantum system, while the k-th element of the corresponding velocity field ΛS (qS , α) is specified according to the elements of the reduced density matrix ΛSk (qS , α)

~ Im = (2π)N mk

(P

0 lS ,lS

σlS0 ,lS (α)ηl∗S (qS ) ∂q∂k ηlS0 (qS ) ρR eq (qS , α)

)

,

(18)

with k ∈ S. We stress that the diffusion matrix does not have components along phase variables α even if they are included amongst the relevant variables xR , since these variables follow a deterministic and independent evolution controlled by the constant velocity field with components {Ej /~}. The underlying reason is that the Schrödinger dynamics, and so 14

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

the quantum phase dynamics, is not affected by the Bohm coordinates even in the stochastic representation that excludes the environment coordinates QB . Consequently, the FokkerPlanck equation (17) enforces a deterministic dynamics for the phases α and a stochastic evolution for the Bohm coordinates of interest qS . Due to the presence of both the streaminglike and the collisional-like term, eq. (17) has the structure of the generic Fokker-Planck equation and it will be denoted by us as the Fokker-Planck-Bohm equation to emphasize its capability to describe the dynamics of the Bohm coordinates in the case of an open quantum system. Indeed, the physical meaning of the Fokker-Planck-Bohm equation is that of i) a purely deterministic dynamics for the phases α driven by the velocity field with components {Ej /~} and ii) a more complex dynamics of the open system coordinates because of the action of both the deterministic drift due to the velocity field ΛS (qS , α) and the collisionallike term that controls their fluctuations.

2.5

Fokker-Planck-Bohm equation driven by the reduced density matrix

The Fokker-Planck-Bohm equation (17) must be considered as the main result of our analysis, since it describes the evolution of Bohm coordinates QS by fully taking into account the Schrödinger dynamics through the phases. This is certainly a valuable feature for accurately reproducing the dynamics of the open quantum system coordinates QS (t), but it represents, at the same time, a significant drawback: the full evolution of the wave function of the whole isolated system, through its phases, must be exactly known even if our interest concerns only the open system. As a matter of fact, this is very demanding from the computational point of view in the case of large isolated systems, since it requires the full diagonalization of the system Hamiltonian. Fortunately, this requirement is not strictly necessary because of the deterministic dynamics of the phases in the Fokker-Planck-Bohm equation. Let us consider the case of the system for a given pure stare as specified by the set of populations and the set A(t) of 15

ACS Paragon Plus Environment

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

Page 16 of 42

the linearly time dependent phases, Aj (t) = Aj (0) + Ej t/~. This corresponds to a time dependent distribution ρR (qS , α, t) on the open system coordinates qS , but with assigned phases α according to a Dirac delta function δ (α − A(t)). Formally such a distribution is recovered as the solution of the Fokker-Planck-Bohm equation (17) for the initial condition ρR (qS , α, 0) = δ (α − A(0)) ρS (qS , 0) with a generic initial distribution ρS (qS , 0) on the open system coordinates qS . The corresponding solution of the Fokker-Planck-Bohm equation is specified as ρR (qS , α, t) = δ (α − A(t)) ρS (qS , t),

(19)

with the following Fokker-Planck equation for the marginal distribution ρS (qS , t) on the open system coordinates, h ∂ S ρ (qS , t) = − ∇S · ΛSσˆ (t) (qS )+ ∂t  −1 i ρS (qS , t) − ∇S · D ρSσˆ (t) (qS )∇S ρSσˆ (t) (qS )

(20)

ˆ S ρS (qS , t), = −Γ σ ˆ (t)

ˆ S acts on the qS coordinates only. An important inwhere the corresponding operator Γ σ ˆ (t) gredient of this time evolution operator is the distribution ρSσˆ (t) (qS ) whose time dependence is parametrically determined by the reduced density matrix σ ˆ (t) = σ ˆ (α)|α=A(t) , and which is obtained from the equilibrium distribution ρR eq (qS , α) by inserting the actual values A(t) of the phases ρSσˆ (t) (qS ) = (2π)N ρR eq (qS , α) α=A(t) X σlS0 ,lS (t)ηl∗S (qS )ηlS0 (qS ). =

(21)

0 lS ,lS

The factor (2π)N in the above equation has been inserted in order to recover a probability density which is normalized by integration on qS coordinates only, i.e., without integration on the phases that was required for the original equilibrium distribution ρR eq (qS , α). It is evident that ρSσˆ (t) (qS ) plays the same role of the equilibrium distribution ρR eq (qS , α) of the Fokker-Planck-Bohm equation (17). The other new ingredient of the eq. (20) is the velocity field ΛSσˆ (t) (qS ) parametrically dependent on the reduced density matrix σ ˆ (t), which again 16

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

derives from the replacement of the phases α with their actual value A(t) for the given pure state in the velocity field ΛS (qS , α) eq. (18) for the Fokker-Planck-Bohm equation ΛSσˆ (t),k (qS ) = ΛSk (qS , α) α=A(t) ). (P ∂ ∗ 0 (qS ) η 0 σl0 ,l (t)ηl (qS ) l ~ lS ,lS ∂q S S S S k = Im mk ρSσ(t) (qS )

(22)

In conclusion, by specializing the general Fokker-Planck-Bohm equation (17) for a given pure state (wave function), we get the time evolution equation (20) for the distribution ρS (qS , t) ˆS on the coordinates qS of the open quantum system, with a time evolution operator Γ σ ˆ (t) which is parametrically dependent on the reduced density matrix σ ˆ (t). Correspondingly, the time dependent reduced density matrix σ ˆ (t) brings all the required information about the quantum state (wave function) that is relevant for the dynamics of the coordinates qS of the open system components. For this reason, we can designate eq. (20) as the Fokker-Planck-Bohm equation driven by the reduced density matrix, even if we omit such a further specification when it is not strictly necessary. Furthermore, we stress that the reduced density matrix plays the role of a pilot field for the coordinates of the open quantum system, very much like the role of the wave function in the deterministic Bohm equation (1). We would like also to mention that the use of the solutions of the Fokker-Planck-Bohm equation (17) for the particular initial condition ρR (qS , α, 0) = δ (α − A(t)) ρS (qS , 0) leading to eq. (20), should not be considered as a limitation of our analysis. Indeed, since any initial density distribution ρR (qS , α, 0) can R be formally written as ρR (qS , α, 0) = dA(0) δ(α − A(0))ρR (qS , A(0), 0), also any solution of eq. (17) can be specified by integration over the initial phases A(0) of the solutions of the

Fokker-Planck-Bohm equation (20) driven by the reduced density matrix. The distribution ρSσˆ (t) (qS ) establishes a clear correspondence between the stochastic description supplied by the Fokker-Planck-Bohm equation and the predictions of standard Quantum Mechanics. Indeed, the same distribution ρSσˆ (t) (qS ) of eq. (21) is obtained from the wave function corresponding to the pure state of eq. (2) by integration on the environmental

17

ACS Paragon Plus Environment

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

degrees of freedom, ρSσˆ (t) (qS )

=

Z

dqB |Ψ(q, t)|2q=(qS ,qB ) .

Page 18 of 42

(23)

In order to emphasize this correspondence between ρSσˆ (t) (qS ) and |Ψ(q, t)|2 after integration on the environment coordinates, we refer to ρSσˆ (t) (qS ) as the quantum marginal distribution of the Fokker-Planck-Bohm equation. It ensures that the average value of any observable of the type b(qS ) according to the quantum marginal distribution ρSσˆ (t) (qS ), Eσˆ (t) [b] = R dqS b(qS )ρSσˆ (t) (qS ), is formally equivalent to the corresponding expectation value hb(qS )it =

hΨ(t)|b(ˆ qS )|Ψ(t)i evaluated according to the pure state |Ψ(t)i.

In this framework, one can also consider the particular situation of an open quantum system in equilibrium with the environment, when the reduced density matrix σ ˆ (t) can be effectively substituted by its time average σ ˆ , since its time evolution is characterized by ˆ . 43,44 The corresponding marginal distribution ρSσˆ (qS ) is time negligible fluctuations around σ independent and it becomes the equilibrium distribution for the open system coordinates qS : ρSeq (qS ) =

X

σ lS0 ,lS ηl∗S (qS )ηlS0 (qS ).

(24)

0 lS ,lS

Furthermore, as already discussed in detail in Ref. 38, a vanishing velocity field ΛσSˆ (qS ) is recovered in correspondence of a Fokker-Planck equation without the streaming-like contribution. Thus, the Fokker-Planck-Bohm equation (20) is reduced to a Smoluchowski type equation

h −1 i S ∂ S ρ (qS , t) = ∇S · D ρSeq (qS )∇S ρSeq (qS ) ρ (qS , t), ∂t

(25)

that is the Smoluchowski-Bohm equation already examined in Ref. 38. In this way, one can shows that the description supplied by the Fokker-Planck-Bohm equation is equivalent to that of the Smoluchowski-Bohm equation when the condition of validity of the later is fulfilled: whenever the time dependence of σ ˆ (t) is characterized by negligible fluctuations around its time average σ ˆ describing the open quantum system in equilibrium with the environment.

18

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

Finally, we report the Langevin equation associated to eq. (20), d QS (t) = ΛSσˆ (t) (QS (t)) + dt −1  D ∇S · ρSσˆ (t) (QS (t)) + (2D)1/2 ξ(t), + ρSσˆ (t) (QS (t))

(26)

that describes a fluctuating stochastic trajectory QS (t) driven by the reduced density matrix.

3

Model system and calculation procedures

In order to compare the deterministic and the stochastic descriptions of the open quantum system, we consider a model system composed of n interacting harmonic oscillators resembling n interacting vibrational degrees of freedom. The set q = (q1 , q2 , . . . , qn ) includes the coordinates of all the oscillators. Because of the similarities between this model system and that analyzed in Refs. 38 and 40, we just summarize its main features. The Hamiltonian of the whole system is specified as ˆ =H ˆ (0) + λVˆ = H

n X

ˆ (0) + λVˆ , H k

(27)

k=1

where Vˆ represents the interactions between the oscillators to be modeled as a random matrix ˆ (0) is the Hamiltonian operator of with the parameter λ measuring their magnitude, and H k the k-th harmonic oscillator, 2 2 ˆ (0) = − ~ ∂ + 1 mk ω 2 q 2 , H k k k 2mk ∂qk2 2

(28)

with ωk (mk ) being its resonance angular frequency (mass). We denote the eigenvalues and ˆ (0) as (0) = (lk + 1/2)~ωk the eigenfunctions of the each oscillator Hamiltonian operator H k lk and ζlk (qk ), respectively. One oscillator (the first one by hypothesis) is taken as the open quantum system of interest S, whereas the other oscillators are considered as the environment B. Correspondingly, the eigenfunctions ζl1 (q1 ) can be employed as basis set for the Hilbert space HS for the open system, ηlS (qS ) = ζl1 (q1 )|l1 =lS ,q1 =qS . The angular frequency and Bohm coordinate of the first oscillator are tagged by ωS and QS (t), respectively. In this case, the open system has 19

ACS Paragon Plus Environment

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

Page 20 of 42

one only degree of freedom whose stochastic dynamics is described by a one-dimensional Fokker-Planck-Bohm equation including a scalar diffusion coefficient D. We employ numerical methods for computing the deterministic and the stochastic evolution of the quantum system. In particular, for solving the Schrödinger equation, we ˆ (0) eigenstates consider the matrix representation of the Hamiltonian on the basis set of H N P (0) (0) |li = nk=1 |ζlk i, where l = (l1 , l2 , . . . , ln ). A diagonal matrix with elements El = nk=1 lk

ˆ (0) , whereas the operator Vˆ has been modelled with a is recovered from the Hamiltonian H random matrix 49,50 whose elements are gaussian distributed according to the following statistical constraints: E [Vl,l0 ] = 0, E [(Vl,l0 )2 ] = 1/(2 − δl,l0 ). Then, the numerical diagonalization ˆ supplies its eigenvalues {Ej } and its eigenvectors {|φj i} that allow one to specify the of H general solution of the Schrödinger equation according to eq. (2), once the initial pure state |Ψ(0)i is established. Let us now consider the problem of selecting the initial wave function in condition of non-equilibrium. We suppose the open quantum system is initially unentangled with the environment and it starts to interact with the other oscillators at t = 0. Correspondingly, the initial wave function |Ψ(0)i is factorized into an open system wave function |ψ S i and an environment wave function |ψ B i, |Ψ(0)i = |ψ S i ⊗ |ψ B i.

(29)

In order to observe the dephasing and the relaxation between the lower energy states of the open system, we choose the wave function |ψ S i as a linear combination of its ground state and its first excited state. The environment wave function |ψ B i is conveniently specified like in eq. (2), as the following linear combination of the eigenstates of the environment P ˆ (0) Hamiltonian nk=2 H k |ψ B i =

X q −ıAB PlB2 ,...,ln e l2 ,...,ln |ζl2 . . . , ζln i,

(30)

l2 ,...,ln

where {PlB2 ,...,ln } and {AB l2 ,...,ln } are the corresponding sets of normalized populations and phases, respectively. With reference to an isolated environment in thermal equilibrium, the 20

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

wave function |ψ B i is selected at random from an uniform distribution on the Bloch sphere P (0) (0) for the active space spanned by the environment eigenstates with energy El2 ,...,ln = nk=2 lk

smaller than a cutoff energy Emax . This correspond to a realization of the so-called Random Pure State Ensemble which has been introduced to characterize the statistical properties of randomly selected quantum states, as well as their thermal properties controlled by cutoff energy Emax . 41,43 In particular the populations {PlB2 ,...,ln } are conveniently selected by employing the algorithm reported in Ref. 51, while the phases {AB l2 ,...,ln } are homogeneously sampled at random in their domain [0 − 2π). For a sufficiently large number of oscillators, the environment behaves like a thermal bath which drives to equilibrium the open system. It should be stressed that with this procedure, the initial wave function |Ψ(0)i is specified ˆ (0) eigenstates |li having (Emax + 3~ωS /2) as upper energy. as a linear combination of H Once the initial pure state |Ψ(0)i is chosen, it is convenient to expand it on the basis ˆ eigenstates |φj i in order to recover a simple representation of of the full Hamiltonian H the dynamics according to eq. (2) in terms of the set P of constant population and of the set A(t) of linear time dependent phases. Correspondingly, also the instataneuos reduced density matrix σ ˆ (t) for the open quantum system is determined according to eq. (13) by ˆ eigenstates |φj i and inserting the phases A(t). This procedure requires the Hamiltonian H eigenvalues Ej , which should be obtained by diagonalization of the matrix representation of ˆ on the oscillator basis |li. The numerical diagonalization, however, has to be performed H on a truncated form of this matrix in correspondence of a finite dimensional Hilbert space (0)

Htr (the “truncated” Hilbert space) generated by basis elements |li with energies El smaller nL o (0) that the truncation cutoff Etr : Htr := |li with E < E . Because of the use of a tr l l

finite basis set, the numerical diagonalization supplies only approximate eigenstates |φj i of the Hamiltonian operator. However, as long as the interaction potential is a perturbation ˆ (0) , the diagonalization of the full Hamiltonian is influenced mainly by the compared to H (0)

coupling amongst basis elements |li with similar energies El . Thus, the basis set truncation has significant effects only for the eigenstates |φj i with eigenenergies close to truncation 21

ACS Paragon Plus Environment

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

energy cutoff Etr . On the other hand, the initial wave function |Ψ(0)i includes |li components (0)

with El

energies up to (Emax + 3~ωS /2). Therefore, in order to deal with a time dependent

wave function that is not affected by the Hamiltonian matrix truncation, one should employ a truncation cutoff Etr sufficentily larger than the (Emax + 3~ωS /2), since in this way the eigenstates affected by numerical truncation do not contribute significantly to |Ψ(0)i. Once the time dependent pure state |Ψ(t)i is determined, the deterministic Bohm trajectory of the ensemble of oscillators is computed by solving the first of eq. (1) according to the Runge-Kutta method 52 at the 4-th order. The stochastic Bohm trajectory of the first oscillator, i.e., the open quantum system, is calculated from the Langevin equation (26) by employing the Euler method, 52 which it is computationally less expensive than the RungeKutta method. Even if the Euler method is unstable in the case of deterministic equations, it results to be an efficient algorithm for stochastic equations because of the white noise term, whose values are sampled according to a gaussian distribution with variance set to the reciprocal of the integration time step.

4

Deterministic dynamics vs stochastic dynamics

In this section, we compare the predictions of our stochastic model with those of the deterministic description according to the Schrödinger-Bohm dynamics. To this aim, we examine first of all the Schrödinger dynamics alone and in particular the time evolution of the reduced density matrix for the open system. Secondly, we investigate the deterministic evolution and the statistical properties of the Bohm coordinate QS (t). Thirdly, we analysis the stochastic dynamics of the first oscillator as supplied by the Fokker-Planck-Bohm equation. In our calculation, we assume that all the oscillators are characterized by the same mass, mk = m ∀k, but by different angular frequencies. In particular, the same angular frequency ω is assigned to the oscillators composing the environment, ωk = ω for k = 2, 3, . . . , n, while the angular frequency of the open system is taken three times larger ωS = 3ω. In this way a

22

ACS Paragon Plus Environment

Page 22 of 42

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

The Journal of Physical Chemistry

sufficiently dense environment levels are available to exchange energy with the open system. The coupling parameter in eq. (27) is set to the value λ = 0.001ω~ to ensure that the interaction plays a perturbative role in coupling the different oscillators. Furthermore, we employ the value of Etr = 12.5ω~ for the truncation cutoff corresponding to a dimension of 1560 for the computational Hilbert space, and the value Emax = 7ω~ for the energy cutoff of the environment states. In this way, the initial wave function |Ψ(0)i is described by a linear ˆ (0) eigenstates |li that have 11.5ω~ as the largest energy to be compared combination of H with the truncation cutoff of 12.5ω~. In the following, we shall employ 2π/ω as time unit p and ~/mωS as the coordinate unit for the first oscillator representing the open system S. One might wonder what would happen by employing different values of the model system

parameters, e.g., the oscillator frequencies, the cutoff energy and the coupling parameter λ. As matter of fact, this would determine a change on the physical properties and the time evolution with respect to the behaviour examined in the next paragraphs. On the other hand, the main results discussed in this work regarding the correspondence between the deterministic and the stochastic description are independent of this choice as long as the model system still represents an open system interacting with a larger environment.

4.1

Schrödinger dynamics

We specifically investigate the case of an initial wave function |ψ S i with equally populated ground and first excited state 1 |ψ S i = √ (|ζ0 i + |ζ1 i), 2

(31)

for the oscillator representing the open quantum system. Then, from the diagonalization of the full Hamiltonian, we determine the time dependent reduced density matrix σ ˆ (t) and its time average σ ˆ . On the basis set of the oscillator eigenstates, |ηlS i = |ζlS i, the matrix representation of σ ˆ results to be almost diagonal, |σ lS ,lS0 |/|σ lS ,lS − σ lS0 ,lS0 | < 1/1000 ∀lS 6= lS0 . In Table 1 the time-averaged diagonal values σ lS ,lS are reported in comparison with the corresponding initial values σlS ,lS (0). The averaged reduced density matrix describes the 23

ACS Paragon Plus Environment

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

Page 24 of 42

Table 1: Diagonal elements of the initial σlS ,lS (0) and of the equilibrium σ lS ,lS reduced density matrix. lS

σlS ,lS (0)

σ lS ,lS

0 1 2

0.5 0.5 0

0.893 0.102 4.15 10−3

equilibrium state, that is the stationary state produced by the equilibration 53 with the set of oscillators representing the environment. In this case the equilibrations leads to deexcitation of the open quantum system with respect to the initial condition, with an increase of the population σ0,0 of the ground state and a corresponding decrease of that σ1,1 of the first excited state, while the population σ2,2 of the following state remains close to zero. Because of the particular conditions chosen for the model system, the Schrödinger dynamics determines a energy flux from the open quantum system to the environment. Given the averaged density matrix σ ˆ , the equilibrium distribution ρSeq (qS ) on the Bohm coordinate of the oscillator is evaluated according to eq. (24). One recovers the confined distribution reported in Fig. 1, which is well reproduced by a normal (gaussian) distribution, 2

2

e−qS /2σqS NσqS (qS ) = √ , 2πσqS

with an optimized width σqS = 0.789

p

(32)

~/mωS . It should be noted that the width is larger

than that of the normal distribution for the ground state of the oscillator, because of the contribution weighted by σ 1,1 of the first excited state. The relaxation dynamics for the populations of the open quantum system is evident by examining the time dependence of the σ0,0 (t) and σ1,1 (t) which are displayed in Fig. 2. Starting from the same value of 0.5, within a time window of 15 (2π/ω) they get very close to the equilibrium value σ 0,0 and σ 1,1 which are reported in Fig. 2 as dashed horizontal line. In the same time window, the corresponding off-diagonal element σ0,1 (t) oscillates, as one can see by examining the time dependence of its real part which is also displayed in Fig. 2. A persistent oscillation with frequency ωS would represent the coherence dynamics of the 24

ACS Paragon Plus Environment

Page 25 of 42

ρSeq (qS ) NσqS (qS )

0.5 0.4 0.3 0.2

p

mωS /¯h, NσqS (qS )/

p

mωS /¯h

0.6

ρSeq (qS )/

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

The Journal of Physical Chemistry

0.1 0.0

−4

−3

−2

1 −1 p 0 ¯ /mωS qS / h

2

3

4

Figure 1: Equilibrium density distribution ρSeq (qS ) of the qS coordinate (solid green line) and p its gaussian approximation NσqS (qS ) with variance σqS = 0.789 ~/mω (dashed black line). isolated oscillator. Instead, the interactions with the environment produces the decoherence, that is the decay of the amplitude of this oscillation as clearly shown in Fig. 2. We mention that the processes of decoherence can be detected also by looking at the time dependence of the expectation value hqS it of the oscillator coordinate qS , because hqS it ∝ Re {σ0,1 (t)} if only the first two states are populated. In conclusion, by examining the Schrödinger dynamics, one can recognize two different regimes. In the first, for times 0 < t < 15 (2π/ω), the oscillator representing the open quantum system clearly displays relaxation of the populations and decoherence for the offdiagonal elements of the reduced density matrix. Afterwards, the oscillator is at equilibrium with the environment, and one can detect only small deviations from the averaged density matrix, in particular the small amplitude oscillations of its off-diagonal elements. The finite size of the environment employed in our calculations is responsible for these deviations from the time average, 43,44 which would be reduced by increasing the number of oscillators for the 25

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1.0 0.8 0.6 0.4

 Re σlS′ ,lS (t)

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

Page 26 of 42

0.2 0.0

−0.2 −0.4

σ0,0(t)

−0.6

Re{σ0,1 (t)}

−0.8

σ1,1(t)

0

5

10

15 ωt/2π

20

25

30

Figure 2: Time evolution of the reduced density matrix elements σ0,0 (red line), σ1,1 (blue line) and of the real part of σ0,1 (green line). The dashed lines represent the corresponding equilibrium values. environment.

4.2

Bohm dynamics

The deterministic evolution of the Bohm coordinates according to eq. (1) has been computed i) under the assumptions that the initial positions of the oscillators are at the bottom of the harmonic potential and ii) by adopting a time step of ∆t = 1.6 · 10−3 2π/ω. We have verified that the trajectories do not change significantly within the correlation time for a further decreases of the time step ∆t of integration. The Bohm trajectory for the open system coordinate QS (t) is displayed in Fig. 3 and one can observe that also the Bohm coordinate evolution is characterized by two different dynamics regimes like the reduce density matrix. Initially the trajectory behaviour is approximately an oscillation at the resonance frequency ωS that becomes at longer times a pure fluctuating dynamics in the presence of a

26

ACS Paragon Plus Environment

Page 27 of 42

p QS (0)/ h ¯ /mωS = 0 p ¯ /mωS = 0.75 QS (0)/ h

3 2 p QS (t)/ h ¯ /mωS

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

The Journal of Physical Chemistry

1 0

−1 −2 −3

0

5

10

15 ωt/2π

20

25

30

Figure 3: Deterministic time evolution of the Bohm coordinate of the first oscillator composing the model system for two different initial coordinates QS (0). confining potential. By comparing Fig. 2 and Fig. 3, it can be verified that the oscillatinglike profile of the Bohm trajectory QS (t) corresponds to the out of equilibrium dynamics of the reduced density matrix σ ˆ (t). On the other hand, pure fluctuating evolution emerges when the Schrödinger dynamics has reached the equilibrium. In this way the Bohm trajectory features are indicative of the corresponding reduced density matrix dynamics. When QS (t) oscillates the reduced density matrix is evolving toward the equilibrium, whereas a fluctuating evolution of QS (t) is symptomatic of the fact that the reduced density matrix has attained its equilibrium values. This behaviour is independent of the initial Bohm coordinate of the oscillator, as one can verify in Fig. 3 where we have also reported a second p Bohm trajectory QS (t) for a different initial condition, i.e., QS (0) = 0.75 ~/mωS . Since the oscillating behaviour characterizes only the initial transient of the trajectory, it does not

influence the statistical properties of the full trajectory, such as the auto-correlation function (and its spectral density) that will be examined later on. 27

ACS Paragon Plus Environment

The Journal of Physical Chemistry

S By sampling a long enough trajectory QS (t), one can evaluate the distribution weq (qS )

on the Bohm coordinate which can be employed to calculate the time average b(QS ) of any observable specified as a function b(QS ) of the Bohm coordinate, Z 1 T dt b (QS (t)) b(QS ) := lim T →+∞ T 0 Z S = dqS b(qS )weq (qS ).

(33)

S (qS ) is indepenIt should be stressed that, because of the long time limit, the distribution weq

dent of the initial transient and it represents the probability density on the Bohm coordinate S (qS ) is displayed in Fig. 4 in comparison to in equilibrium conditions. The distribution weq

the equilibrium distribution ρSeq (qS ). The equivalence of the two distributions is evident, de-

S weq (qS )

ρSeq (qS )

0.5 0.4 0.3 0.2

p

mωS /¯h, ρSeq (qS )/

p

mωS /¯h

0.6

S weq (qS )/

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

Page 28 of 42

0.1 0.0

−4

−3

−2

1 −1 p 0 ¯ /mωS qS / h

2

3

4

S Figure 4: Distribution weq (qS ) (red line) of the Bohm coordinate corresponding to the first oscillator calculated from the sampling of its deterministic trajectory QS (t). For comparison the equilibrium distribution ρSeq (qS ) (green line) derived from the reduced density matrix is also reported.

S spite the roughness of the profile of weq (qS ) due to the finite time window T of the calculated

trajectory QS (t) (in this case T = 6000 2π/ω). This represents an important result: the two 28

ACS Paragon Plus Environment

Page 29 of 42

distributions derive from different kinds of information, ρSeq (qS ) from the wave function and S (qS ) from the Bohm trajectory, so that their equivalence should not be taken for granted. weq

This issue has been examined in detail by us in a previous work 7 and it can been explained on the basis of the ergodic properties of the Bohm dynamics. Let us determine the time scale of the fluctuations which are displayed by the evolution of the Bohm coordinate QS (t). To this aim, we employ its auto-correlation functions, GQS (τ ) :=

∆QS (t + τ )∆QS (t) (∆QS (t))2

,

(34)

where the overline denotes the time average like in eq. (33) and ∆QS (t) = QS (t)−QS (t). For a given time delay τ , it is evaluated by averaging directly ∆QS (t+τ )∆QS (t) for 0 ≤ t ≤ T −τ along the trajectory QS (t) calculated for a large enough time window T . Again, as long as

1.0 deterministic stochastic Ornstein-Uhlenbeck

0.8 0.6 GQS (τ )

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

The Journal of Physical Chemistry

1.00

0.4 0.2

0.99 0.000

0.025

0.0 0

5

10

15 ωτ /2π

20

25

30

Figure 5: Correlation function GQS (τ ) of the coordinate QS (t) of the first oscillator from the deterministic trajectory (red line); from the stochastic trajectory corresponding to the Fokker-Planck-Bohm equation (black line); from the Ornstein-Uhlenbeck process (dashed black line). The inset magnifies the initial decay in order to recognize the different time parity of the deterministic and stochastic time evolution. we employ a large time window T , the resulting function GQS (τ ) is independent of the initial 29

ACS Paragon Plus Environment

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

transient and it describes the correlation of the QS (t) fluctuations in the equilibrium regime. In Fig. 5, we display with the red profile the calculated correlation function together with the exponential decay GQS (τ ) ∝ e−τ /τC obtained from the best fit (dotted line in Fig. 5), which allows us to estimate the correlation time of QS fluctuations as τC = 1.87 2π/ω. It is worthwhile to spend few more words for what concerns the effects of the coupling parameter λ in eq. (27) on the correlation time τC . Indeed, we have found that by either decreasing or increasing the strength parameter λ (within the boundaries of perturbative interactions), the correlation time is not modified. This suggests that the correlation time of the Bohm coordinates is mainly controlled by the entanglement between the open system and the environment instead of the magnitude of the interaction. At the actual stage, this trend has to be considered as a numerical observation that calls for further investigations.

4.3

Stochastic dynamics

We now analyze the Bohm dynamics by employing the stochastic representation according to the Fokker-Planck-Bohm equation (20) or the corresponding Langevin type equation (26). To this aim, three fundamental ingredients, i.e., the quantum marginal distribution ρSσˆ (t) (qS ), the velocity field ΛSσˆ (t) (qS ) and the diffusion coefficient D, have to be determined. As regards the quantum marginal distribution ρSσˆ (t) (qS ) and the velocity field ΛSσˆ (t) (qS ), they can be exactly computed through eq. (21) and eq. (22), respectively, once the time dependent reduced density matrix σ ˆ (t) is determined. As already discussed, for our model system σ ˆ (t) is exactly calculated according to eq. (13) and the time evolution of its elements is displayed in Fig. 2. Let us now examine the issue of determining the diffusion coefficient D. Similarly to what is often done in the framework of classical stochastic methods, we have assumed that the diffusion coefficient does not depend on the relevant variables. Then the same diffusion coefficient can be used for the initial transient and the equilibrium regime described by the Smoluchowski-Bohm equation (25). This allows us to determine the value of D by 30

ACS Paragon Plus Environment

Page 30 of 42

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

The Journal of Physical Chemistry

comparing the autocorrelation function GQS (τ ) calculated with the deterministic dynamics and displayed in Fig. 5 with the prediction of the Smoluchowski-Bohm equation (25) for the same autocorrelation function. For the former case we have obtained the correlation time τC = 1.87 2π/ω. On the other hand, the mono-exponential decay of the OrnsteinUhlenbeck process 18 with τC = σq2S /D is derived for the same autocorrelation function from the Smoluchowski-Bohm equation 38 if the normal (gaussian) distribution of eq. (32) is used for the equilibrium probability density ρSeq (qS ). Then, given the equilibrium variance p σqS = 0.789 ~/mωS of the oscillator coordinate that ensures the best fit of ρSeq (qS ) with the

normal distribution as displayed in Fig. 1, we obtain the estimate D = 0.0529 ~/m for the

diffusion coefficient, which has been used for the following stochastic calculations. The comparison between the stochastic trajectory QS (t) resulting for the Langevin equation (26) and its deterministic counterpart starting at the same point QS (0) = 0 and which has been already examined in Fig. 3, is made in Fig. 6. For the numerical solution of the Langevin equation, we have employed a time step of ∆t = 8.0 · 10−3 2π/ω. Clearly, the two trajectories share the same features: the initial dynamics is roughly an oscillation at the resonance frequency of the oscillator ωS , which progressively becomes a purely fluctuating evolution when the equilibrium regime is attained. The agreement of the resulting behaviour does not allow one to recognize between the two trajectories which is the stochastic or the deterministic one by a direct inspection of Fig. 6. Similar differences can be found between the two deterministic trajectories of Fig. 3 and, therefore, one could interpret the stochastic trajectory of Fig. 6 as another deterministic trajectory deriving from a different initial condition. We would like to stress that the initial oscillating behaviour can not be recovered from the Smoluchowski-Bohm equation (25), since it is an effect of the streaming-like contribution due to the velocity field ΛSσ(t) (qS ) in the Fokker-Planck-Bohm representation eq. (20). The comparison in Fig. 5, between the deterministic and the stochastic correlations functions, highlights the effectiveness of the projection procedure leading to the Fokker-PlanckBohm equation and of the estimate of the diffusion coefficient as well. Indeed both the cor-

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry

deterministic stochastic

3 2 1

p QS (t)/ h ¯ /mω

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

Page 32 of 42

0

−1 −2 −3

0

5

10

15 ωt/2π

20

25

30

Figure 6: Stochastic evolution (black line) of the coordinate of the first oscillator as obtained from the Langevin equation. For comparison, the deterministic trajectory of Fig. 3 with the same initial condition (QS (0) = 0) is also reported (red line). relation functions emerging from the Langevin dynamics and from the Ornstein-Uhlenbeck process, i.e., the exponential decay, well reproduce the profile of the deterministic correlation function. However, one should take into account that projection procedures generates a stochastic representation which is necessarily an approximation of the deterministic dynamics. The aim of the stochastic representation is that of reproducing the main features of the deterministic motion without a truly quantitative agreement and, therefore, deviations should be expected. For instance, differences between the two kinds of autocorrelation functions emerge at short times. The deterministic dynamics, being time reversible, leads to a correlation function with a even time parity so that dGQS (τ )/dτ |τ =0 = 0. On the other hand, in stochastic models having a semi-group evolution describing dissipation phenomena, this reversibility is absent and the corresponding correlation functions are no more even functions of the time, i.e., dGQS (τ )/dτ |τ =0 6= 0. Such a difference is evident in the inset in Fig. 5 32

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

which magnifies the correlation function in a time domain much shorter than the correlation time. However, the different time parity does not have a significant effect on the comparison between the correlation functions in a time window of the order of the correlation time τC . Finally, we would like to show how the stochastic representation can be employed to calculate time dependent quantum expectation values by means of a swarm of stochastic trajectories. As already discussed, for any observable specified as a function b(qS ) of the coordinates, the expectation value hb(qS )it = hΨ(t)|b(ˆ qS )|Ψ(t)i evaluated according to the Schrödinger dynamics of the the pure state |Ψ(t)i is formally equivalent to the average R Eσˆ (t) [b(qS )] = dqS b(qS )ρSσˆ (t) (qS ) according to the quantum marginal distribution ρSσˆ (t) (qS ) of the Fokker-Planck-Bohm equation. This can be numerically verified by averaging the

observable on a swarm of Langevin trajectories whose initial conditions QS (0) are sampled according to quantum marginal distribution ρSσˆ (t) (qS )|t=0 . We report in Fig. 7 the comparison between these two different kinds of calculations in the case of the open system coordinate, i.e., b(qS ) = qS . The calculation have been performed with 10000 trajectories. The two profiles are so perfectly overlapped that, in order to recognized them, we have employed in Fig. 7 two separates plots. This demonstrates the capability of the Fokker-Planck-Bohm equation to reproduce correctly both the dynamical process characterizing the coordinate evolution of the open quantum system and the quantum predictions according to the Schrödinger dynamics.

5

Conclusions

In this work, we have developed a stochastic formulation of the quantum molecular trajectory in order i) to represent accurately the molecular motions of an open system within a quantum framework and ii) to avoid the high computational cost of the deterministic dynamics. Indeed, the quantum molecular trajectory supplies a complete and deterministic description of the molecular motions, but it requires the solution of both the Schrödinger

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry

~/mωS

1

hqS it /

p

0

1

p

~/mωS

−1

Eσˆ(t) [qS ]/

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

Page 34 of 42

0

−1

0

5

10 ωt/2π

15

20

Figure 7: Time evolution of the expectation value for the open quantum system coordinate qS . The first profile, hqS it , has been determined through the exact solution of the Schrödinger equation (geen line), whereas the other, Eσˆ (t) [qS ], has been calculated by averaging qS on a swarm of Langevin trajectories sampled according to ρSσˆ (t) (qS )|t=0 (black line). equation and the Bohm equation for the isolated system. 7 As a matter of fact, this is an unfeasible task if one considers systems of Chemical interest with a large number of degrees of freedom, like in the case of a molecule interacting with the environment. Therefore, one has to use a stochastic representation of selected coordinates of interest, and in this work we have derived the corresponding Fokker-Planck-Bohm equation (20) through projection operator techniques. 34–37 The Fokker-Planck-Bohm equation describes the time evolution of the Bohm coordinates of an open quantum system with the reduced density matrix, determined from the wave function by tracing out the environment, playing the role of a pilot field. This leads to an accurate characterization of the molecular motions whether the open system is in equilibrium with the environment or not. In the former case, the Fokker-Planck-Bohm equation becomes equivalent to the simpler Smoluchowski-Bohm equation that has been introduced and analyzed in detail in Ref. 38. In the case of the out of equilibrium dynamics, 34

ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

the full Fokker-Planck-Bohm equation has to be employed and in this work the accuracy of its predictions has been verified by comparison with the deterministic quantum molecular trajectory for a model system of six interacting harmonic oscillators where one oscillator has been taken as the open system of interest. The intrinsic benefit of this approach is that of characterising the molecular motions according to the stochastic description of few selected degrees of freedom of interest. Because of the generality of the method which is well suited for systems far from equilibrium, it could find useful applications in the molecular interpretation of relevant chemical processes like photoinduced conformational changes (as the photoisomerization of azobenzene) or photoinduced proton transfers, up to consider kinetic processes driven by external field. In order to employ the Fokker-Planck-Bohm equation in a realistic representation of molecular systems, two fundamental ingredients must be determined in advance: the time dependence of the reduced density matrix σ ˆ (t) and the diffusion coefficient D. In our calculations for the model system, we have computed exactly the time evolution of σ ˆ (t), because the limited dimension of the environment allowed us to solve the Schrödinger equation of the isolated system. For an open system of chemical interest interacting with a macroscopic environment this approach is obviously unfeasible. However, one could model the dynamics of the reduced density matrix by employing well known methods such as the Redfield or the Lindblad forms of quantum master equation. 33,54–57 As regards the diffusion coefficient D which determines the time scale of the fluctuating dynamics, in the applications it can not be evaluated according to the deterministic evolution as we have done in this work. Moreover, the quantum counterpart of Stokes-Einstein relation in Classical Mechanics is not available to supply estimates of D on the basis of the physical properties if the systems. Indeed, the clear understanding of the mechanisms determining the value of the diffusion coefficient for the Bohm coordinates is still an open issue and it deserves further investigations. Some information could be obtained by examining model systems, like that considered in this work, in order to capture the trends of the diffusion coefficient

35

ACS Paragon Plus Environment

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

D with respect to the main features of the systems, e.g., the number of components and the magnitude of the interactions with the environment. At the actual stage, it has to be considered as a parameter of the stochastic model, whose values should be tuned on the basis of the predicted behaviours for the the coordinates of the open quantum system.

Acknowledgement The authors acknowledge the support by Unirvesità degli Studi di Padova through 60% grants and the project P-DiSC#08BIRD2016-UNIPD.

References (1) Craig, I. R.; Manolopoulos, D. E. Chemical Reaction Rates From Ring Polymer Molecular Dynamics. J. Chem. Phys. 2005, 122, 084106. (2) Craig, I. R.; Manolopoulos, D. E. A Refined Ring Polymer Molecular Dynamics Theory of Chemical Reaction Rates. J. Chem. Phys. 2005, 123, 034102. (3) Wang, H.; Sun, X.; Miller, W. H. Semiclassical Approximations for the Calculation of Thermal Rate Constants for Chemical Reactions in Complex Molecular Systems. J. Chem. Phys. 1998, 108, 9726–9736. (4) Cao, J.; Voth, G. A. The Formulation of Quantum Statistical Mechanics Based on the Feynman Path Centroid Density. II. Dynamical Properties. J. Chem. Phys. 1994, 100, 5106–5117. (5) Jang, S.; Voth, G. A. A Derivation of Centroid Molecular Dynamics and other Approximate Time Evolution Methods for Path Integral Centroid Variables. J. Chem. Phys. 1999, 111, 2371–2384.

36

ACS Paragon Plus Environment

Page 36 of 42

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

The Journal of Physical Chemistry

(6) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. Chem. Rev. 2016, 116, 7529– 7550, PMID: 27049513. (7) Avanzini, F.; Moro, G. J. Quantum Molecular Trajectory and Its Statistical Properties. J. Phys. Chem. A 2017, 121, 5352–5360, PMID: 28650172. (8) Bohm, D. A Suggested Interpretation of the Quantum Theory in Terms of “Hidden” Variables. I. Phys. Rev. 1952, 85, 166–179. (9) Bohm, D. A Suggested Interpretation of the Quantum Theory in Terms of “Hidden” Variables. II. Phys. Rev. 1952, 85, 180–193. (10) Holland, P. R. The Quantum Theory of Motion; Cambridge University Press, 1995; Cambridge Books Online. (11) Lopreore, C. L.; Wyatt, R. E. Quantum Wave Packet Dynamics with Trajectories. Phys. Rev. Lett. 1999, 82, 5190–5193. (12) Wyatt, R. E. Quantum Dynamics with Trajectories - Introduction to Quantum Hydrodynamics; Springer, 2005. (13) Gindensperger, E.; Meier, C.; Beswick, J. A. Mixing Quantum and Classical Dynamics Using Bohmian Trajectories. J. Chem. Phys. 2000, 113, 9369–9372. (14) Garashchuk, S.; Dell’Angelo, D.; Rassolov, V. A. Dynamics in the Quantum/Classical Limit Based on Selective Use of the Quantum Potential. J. Chem. Phys. 2014, 141, 234107. (15) Figalli, A.; Klein, C.; Markowich, P.; Sparber, C. WKB Analysis of Bohmian Dynamics. Comm. Pure Appl. Math. 2014, 67, 581–620.

37

ACS Paragon Plus Environment

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

(16) Berndl, K.; Dürr, D.; Goldstein, S.; Peruzzi, G.; Zanghì, N. On the Global Existence of Bohmian Mechanics. Commun. Math. Phys. 1995, 173, 647–673. (17) Teufel, S.; Tumulka, R. Simple Proof for Global Existence of Bohmian Trajectories. Commun. Math. Phys. 2005, 258, 349–365. (18) Gardiner, C. W. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences; Springer-Verlag: New York, 1986. (19) Gillespie, D. T. A Rigorous Derivation of the Chemical Master Equation. Physica A 1992, 188, 404 – 425. (20) Gillespie, D. T. Stochastic Simulation of Chemical Kinetics. Annu. Rev. Phys. Chem. 2007, 58, 35–55. (21) Gillespie, D. T.; Hellander, A.; Petzold, L. R. Perspective: Stochastic Algorithms for Chemical Kinetics. J. Chem. Phys. 2013, 138, 170901. (22) Wu, F.; Tian, T.; Rawlings, J. B.; Yin, G. Approximate Method for Stochastic Chemical Kinetics with Two-Time Scales by Chemical Langevin Equations. J. Chem. Phys. 2016, 144, 174112. (23) Prinz, J.-H.; Chodera, J. D.; Noé, F. Spectral Rate Theory for Two-State Kinetics. Phys. Rev. X 2014, 4, 011020. (24) Hirata, F.; Akasaka, K. Structural Fluctuation of Proteins Induced by Thermodynamic Perturbation. J. Chem. Phys. 2015, 142, 044110. (25) Klein, H. C. R.; Schwarz, U. S. Studying Protein Assembly with Reversible Brownian Dynamics of Patchy Particles. J. Chem. Phys. 2014, 140, 184112. (26) Ryabov, Y. E.; Fushman, D. A Model of Interdomain Mobility in a Multidomain Protein. J. Am. Chem. Soc. 2007, 129, 3315–3327. 38

ACS Paragon Plus Environment

Page 38 of 42

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

The Journal of Physical Chemistry

(27) Wong, V.; Case, D. A.; Szabo, A. Influence of the Coupling of Interdomain and Overall Motions on NMR Relaxation. Proc. Natl. Acad. Sci. USA 2009, 106, 11016–11021. (28) Ryabov, Y.; Clore, G. M.; Schwieters, C. D. Coupling Between Internal Dynamics and Rotational Diffusion in the Presence of Exchange Between Discrete Molecular Conformations. J. Chem. Phys. 2012, 136, 034108. (29) Sanz, A.; Martínez-Casado, R.; Peñate-Rodríguez, H.; Rojas-Lorenzo, G.; MiretArtés, S. Dissipative Bohmian Mechanics within the Caldirola–Kanai Framework: A trajectory Analysis of Wave-Packet Dynamics in Viscid Media. Ann. Phys. 2014, 347, 1 – 20. (30) Garashchuk, S.; Dixit, V.; Gu, B.; Mazzuca, J. The Schrödinger Equation with Friction from the Quantum Trajectory Perspective. J. Chem. Phys. 2013, 138, 054107. (31) Polychronakos, A. P.; Tzani, R. Schrödinger Equation for Particle with Friction. Phys. Lett. B 1993, 302, 255 – 260. (32) Kostin, M. D. On the Schrödinger–Langevin Equation. J. Chem. Phys. 1972, 57, 3589– 3591. (33) Breuer, H.; Petruccione, F. The Theory of Open Quantum Systems; OUP Oxford, 2007. (34) Nakajima, S. On Quantum Theory of Transport Phenomena: Steady Diffusion. Prog. Theor. Phys. 1958, 20, 948–959. (35) Zwanzig, R. Ensemble Method in the Theory of Irreversibility. J. Chem. Phys. 1960, 33, 1338–1341. (36) Zwanzig, R. On the Identity of Three Generalized Master Equations. Physica 1964, 30, 1109–1123. (37) Henderson, D. Physical Chemistry: An Advanced Treatise - Mathematical Methods; Elsevier Science, 2012; Vol. XIB. 39

ACS Paragon Plus Environment

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

(38) Avanzini, F.; Moro, G. J. Quantum Stochastic Trajectories: the Smoluchowski-Bohm Equation. Phys. Chem. Chem. Phys. 2018, 20, 165–179. (39) Wyatt, R. E. Quantum Wavepacket Dynamics with Trajectories: Wavefunction Synthesis along Quantum Paths. Chem. Phys. Lett. 1999, 313, 189 – 197. (40) Avanzini, F.; Fresch, B.; Moro, G. J. Pilot-Wave Quantum Theory with a Single Bohm’s Trajectory. Found. Phys. 2016, 46, 575–605. (41) Fresch, B.; Moro, G. J. Emergence of Equilibrium Thermodynamic Properties in Quantum Pure States. I. Theory. J. Chem. Phys. 2010, 133, 034509. (42) Sinai, Y. G. Dynamical Systems II: Ergodic Theory with Applications to Dynamical Systems and Statistical Mechanics; Encyclopaedia of Mathematical Sciences; Springer Berlin Heidelberg, 1989. (43) Fresch, B.; Moro, G. J. Beyond Quantum Microcanonical Statistics. J. Chem. Phys. 2011, 134, 054510. (44) Fresch, B.; Moro, J. G. Typical Response of Quantum Pure States. Eur. Phys. J. B. 2013, 86, 1–13. (45) Efthymiopoulos, C.; Kalapotharakos, C.; Contopoulos, G. Nodal Points and the Transition from Ordered to Chaotic Bohmian Trajectories. J. Phys. A 2007, 40, 12945. (46) Wu, H.; Sprung, D. Quantum Chaos in Terms of Bohm Trajectories. Phys. Lett. A 1999, 261, 150 – 157. (47) de Alcantara Bonfim, O. F.; Florencio, J.; Sá Barreto, F. C. Quantum Chaos in a Double Square Well: An Approach Based on Bohm’s View of Quantum Mechanics. Phys. Rev. E 1998, 58, 6851–6854. (48) Frisk, H. Properties of the Trajectories in Bohmian Mechanics. Phys. Lett. A 1997, 227, 139 – 142. 40

ACS Paragon Plus Environment

Page 40 of 42

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

The Journal of Physical Chemistry

(49) Wigner, E. P. Random Matrices in Physics. SIAM Review 1967, 9, 1–23. (50) Brody, T. A.; Flores, J.; French, J. B.; Mello, P. A.; Pandey, A.; Wong, S. S. M. Random-Matrix Physics: Spectrum and Strength Fluctuations. Rev. Mod. Phys. 1981, 53, 385–479. (51) Fresch, B.; Moro, G. J. Typicality in Ensembles of Quantum States: Monte Carlo Sampling versus Analytical Approximations. J. Phys. Chem. A 2009, 113, 14502– 14513. (52) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes; Cambridge University Press, 2007. (53) García-Pintos, L. P.; Linden, N.; Malabarba, A. S. L.; Short, A. J.; Winter, A. Equilibration Time Scales of Physically Relevant Observables. Phys. Rev. X 2017, 7, 031027. (54) Redfield, A. On The Theory of Relaxation Processes. IBM J. Res. Devel. 1957, 1, 19–31. (55) Redfield, A. In Advan. Magn. Reson; Waugh, J. S., Ed.; Advances in Magnetic and Optical Resonance; Academic Press, 1965; Vol. 1; pp 1 – 32. (56) Lindblad, G. On the Generators of Quantum Dynamical Semigroups. Commun. Math. Phys. 1976, 48, 119–130. (57) Gorini, V.; Kossakowski, A.; Sudarshan, E. C. G. Completely Positive Dynamical Semigroups of N-Level Systems. J. Math. Phys. 1976, 17, 821–825.

41

ACS Paragon Plus Environment

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

6

TOC Graphic

Isolated System

Open System

Interactions with environment lead to decoherence and relaxation processes

QS (t): Quantum Stochastic Trajectories ˆ (t) : Reduced Density Matrix D rS · ⇢Sˆ (t) (QS (t)) d QS (t) = ⇤Sˆ (t) (QS (t)) + + (2D)1/2 ⇣(t) dt ⇢Sˆ (t) (QS (t))

42

ACS Paragon Plus Environment

Page 42 of 42