Quantum Markovian master equation for scattering from surfaces Haifeng Li, Jiushu Shao, Asaf Azuri, Eli Pollak, and Robert Alicki Citation: The Journal of Chemical Physics 140, 014104 (2014); doi: 10.1063/1.4851075 View online: http://dx.doi.org/10.1063/1.4851075 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/1?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Diffusive behavior from a quantum master equation J. Math. Phys. 52, 083303 (2011); 10.1063/1.3614779 The non-Markovian quantum master equation in the collective-mode representation: Application to barrier crossing in the intermediate friction regime J. Chem. Phys. 123, 204111 (2005); 10.1063/1.2121649 Controllability properties for finite dimensional quantum Markovian master equations J. Math. Phys. 44, 2357 (2003); 10.1063/1.1571221 Erratum: “Fourth-order quantum master equation and its Markovian bath limit” [J. Chem. Phys. 116, 2705 (2002)] J. Chem. Phys. 117, 10428 (2002); 10.1063/1.1519534 Fourth-order quantum master equation and its Markovian bath limit J. Chem. Phys. 116, 2705 (2002); 10.1063/1.1445105
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
THE JOURNAL OF CHEMICAL PHYSICS 140, 014104 (2014)
Quantum Markovian master equation for scattering from surfaces Haifeng Li,1 Jiushu Shao,1 Asaf Azuri,2 Eli Pollak,2,a) and Robert Alicki2,b) 1 2
College of Chemistry, Beijing Normal University, Beijing 100875, China Chemical Physics Department, Weizmann Institute of Science, 76100 Rehovoth, Israel
(Received 3 October 2013; accepted 30 November 2013; published online 6 January 2014) We propose a semi-phenomenological Markovian Master equation for describing the quantum dynamics of atom-surface scattering. It embodies the Lindblad-like structure and can describe both damping and pumping of energy between the system and the bath. It preserves positivity and correctly accounts for the vanishing of the interaction of the particle with the surface when the particle is distant from the surface. As a numerical test, we apply it to a model of an Ar atom scattered from a LiF surface, allowing for interaction only in the vertical direction. At low temperatures, we find that the quantum mechanical average energy loss is smaller than the classical energy loss. The numerical results obtained from the space dependent friction master equation are compared with numerical simulations for a discretized bath, using the multi-configurational time dependent Hartree methodology. The agreement between the two simulations is quantitative. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4851075] I. INTRODUCTION
The quantum dynamics of atom surface scattering presents a challenge to theory.1–6 It is relatively easy to solve for the quantum dynamics in the presence of a frozen surface. In this case, one has at most to undertake a three dimensional wavepacket propagation,7, 8 and this is routinely feasible, using, for example, split operator techniques.9 If however, one includes the interaction of the scattered particle with the surface, the problem becomes much more difficult. Even if the interacting surface modes are considered to be harmonic and the interaction is restricted to be linear in the coordinates of the surface modes, the problem is no longer straightforward due to the addition of the thermally distributed continuum of surface modes. In contrast to present day tractable problems10 such as the spin-boson model11, 12 or a system coupled bilinearly to a harmonic bath which may be treated using a variety of methods,13–15 here necessarily, the coupling is only linear in the bath coordinates but it must be nonlinear in the system coordinates. When the particle is far away from the surface the interaction vanishes and it is only felt when the particle is sufficiently close to the surface. It is this nonlinearity which makes a quantum solution for the dynamics difficult. One may well ask—who cares? Typically, if the particle is heavy, then the dynamics may be considered to be classical,16 while if it is sufficiently light, such as He one may use perturbation theory techniques6, 17 due to the very weak interaction of the atom with any surface. In fact though this is not sufficient. Quantum diffraction effects have been measured for heavy atoms and molecules scattered from surfaces.18 In most cases, the resulting angular distributions especially at very low surface temperatures have not received an adequate theoretical analysis. Moreover, it has been shown a) Electronic mail: [email protected] b) Permanent address: Institute of Theoretical Physics and Astrophysics,
University of Gda´nsk, Gda´nsk, Poland.
0021-9606/2014/140(1)/014104/10/$30.00
recently that one may observe quantum effects for heavy particle scattering simply by suitable collimation of the incident beam19, 20 provided that the de Broglie wavelength in the direction perpendicular to the motion of the incident particle is longer than the lattice length of the surface. There are also some good theoretical reasons for attempting to probe deeper into the quantum dynamics of surface scattering. For example, we do not know very much about the difference between the classical and quantum mechanical average energy lost by the particle when colliding with a surface.21 The same is true for the sticking probability to the surface as a function of incident energy or surface temperature.22 The central theme of this paper is to provide a framework for the study of these questions, based on a master equation description of the scattering dynamics. The formulation of a master equation for scattering is not a new problem. An early attempt to develop a formalism for dissipative scattering in terms of Markovian master equations was presented in the context of heavy-ion collisions in Ref. 23. A different approach was developed by Gao and co-workers24, 25 who introduce a master equation which is appropriate in the high-temperature regime as may be inferred by comparison with a generalization of the Caldeira-Leggett quantum Fokker-Planck type of equations. The theory used by Nest26 is of special relevance to the formulation presented in this paper. Using the formalism of super-symmetric quantum mechanics, Nest employed generalized raising and lowering operators for the system Hamiltonian. These operators encode the discrete structure of the Hamiltonian spectrum however they do not include the continuum part, which is inherent to the scattering process. Finally, we mention the class of master equations for a free particle interacting with a homogeneous equilibrium environment.27, 28 These equations satisfy the elementary conditions needed to describe the dissipative dynamics of a continuous spectrum Hamiltonian (see also below).
140, 014104-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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-2
Li et al.
In this paper we formulate a Lindblad like master equation, which introduces generalized Lindblad operators which are space dependent. They vanish when the particle is sufficiently far removed from the surface. In principle, the standard Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) form29, 30 is applicable under the condition that the motion of the system is relatively slow in comparison with the fast dynamics of a bath. The consistent derivation of such equations is available for the cases of well-separated time-scales, i.e., for spatially confined systems with a discrete spectrum of the Hamiltonian and in the weak coupling regime where dissipation is relatively slow in comparison with the system’s Hamiltonian motion. If this is not the case, different Markovian/nonMarkovian regimes are expected which depend on the time and energy scales determined by the initial conditions.31 As already noted, the scattering problems involve initial states (wave packets) composed of energy eigenstates from the continuous part of the system’s Hamiltonian spectrum. Therefore, we are not able to present a rigorous derivation but rather arguments based on extrapolation of the known results for the discrete spectrum case to our model and on the phenomenological properties of the solutions of master equations in the semiclassical regime. Since quantum effects will usually appear at low surface temperatures, it is important to verify that the theory presented is indeed valid in this regime. For this purpose, we compare the zero surface temperature results obtained from the master equation with a full quantum dynamics simulation, using a discretized bath and the Multi-Configurational Time Dependent Hartree (MCTDH) Heidelberg package.10, 32–44 We find excellent agreement, validating the usefulness of the master equation approach. In Sec. II we present the generalized quantum master equation. A numerical example, using parameters which are relevant for the scattering of an Ar atom from a LiF surface is presented in Sec. III. We find, that the average energy lost by the particle to the surface is especially at low temperatures, somewhat lower than the classical. In Sec. IV we describe the adaptation of the MCTDH methodology to the model studied with the master equation and compare the numerical results obtained at T = 0 K. We end with a discussion of the present methodology, its strengths, and weaknesses as well as its further generalization for treating the full three dimensional scattering problem.
II. QUANTUM MARKOVIAN MASTER EQUATIONS FOR DISSIPATIVE SCATTERING A. One-dimensional surface scattering model
In the following we will consider the dynamics of a particle with mass M incident on a smooth surface, interacting with it in the vertical direction only. This model may be readily generalized to include all three degrees of freedom (the vertical z and horizontal x and y directions), however the essential concepts are more readily demonstrated by considering only the vertical motion. The model Hamiltonian (with vertical momentum p) for describing the interaction of the system
J. Chem. Phys. 140, 014104 (2014)
with a harmonic bath takes the form16 ⎡ 2 ⎤ √ N Mc 1 ⎣ 2 p2 j +V (z)+ pj + ωj2 xj − H= V (z) ⎦ . 2M 2 j =1 ωj2 (2.1) The system potential V (z) vanishes for large vertical distances, it has the qualitative form of a Morse potential. The bath Hamiltonian using mass weighted coordinates (xj ) and momenta (pj ) is
1 2 p + ωj2 xj2 , 2 j =1 j N
HB =
(2.2)
while the system-bath interaction term is N √ HI = − MV (z) cj xj .
(2.3)
j =1
The continuum limit for the interaction with the surface is obtained by introducing the spectral density of bath modes: J (ω) =
N 2 π cj δ(ω − ωj ). 2 j =1 ωj
(2.4)
The friction function is then obtained through the Fourier transform of the spectral density 2 ∞ J (ω) γ (t) = . (2.5) dω cos(ωt) π 0 ω Since we want to approximate the dynamics with a Markovian master equation, we must consider the memoryless Ohmic spectral density J (ω) = γ ω
(2.6)
γ (t) = 2γ δ(t)
(2.7)
or conversely
identifying γ as the friction coefficient and δ(x) is the Dirac “delta” function. In the classical limit, with the Ohmic spectral density, the dynamics of the Hamiltonian given in Eq. (2.1) reduces to a space dependent Langevin equation of motion: √ (2.8) M z¨ + V (z) + Mγ [V (z)]2 z˙ = MV (z)F (t), where the (discretized) Gaussian random force is identified as
N pj F (t) = cj xj cos(ωj t) + sin(ωj t) . (2.9) ωj j =1 When averaged over the thermal bath distribution (exp(−βHB ) with β = 1/(kB T)) the random force vanishes, and its thermally averaged autocorrelation function is readily seen to satisfy the classical fluctuation dissipation relation: 1 F (t)F (t ) = γ (t − t ), (2.10) β with 1/β = kB T. Finally, we note that in the limit of weak friction, using classical perturbation theory one readily finds that the leading
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-3
Li et al.
J. Chem. Phys. 140, 014104 (2014)
order term for the average energy loss of the particle after colliding with the surface (E = Ei − Ef , where Ei and Ef are the initial and final energies of the scattered particle) is to first order in the friction ∞ 0 2 dV (zt ) γ ∞ 0 2 Eβ = Mγ dt − dt V zt , dt β −∞ −∞ (2.11) where zt0 denotes the solution of the unperturbed equation of motion M z¨ t0 + V (zt0 ) = 0 and ···β denotes the average with respect to the thermal state of the bath. The first term on the right hand side gives the average energy lost to the bath at T = 0 K, while the second term expresses the fact that as the temperature increases, the energy loss decreases. The bath can now also transfer energy back to the system.
B. LGKS equations for weakly coupled systems
Before writing down the reduced equation of motion for the scattering problem, we review briefly the structure of the Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) equation45, 46 in the weak coupling regime. The reduced equation of motion is constructed such that it preserves the following conditions: (i) The solution of the reduced equation of motion is given by a completely positive, trace preserving one-parameter semigroup. (ii) The diagonal elements of the density matrix in the energy basis evolve independently of the non-diagonal ones. (iii) The thermal canonical distribution is a stationary state. (iv) The detailed balance condition holds. The Schrödinger picture reduced dynamics of the system, in a weak coupling regime, is then given by the following Markovian master equation for the reduced density operator ρˆ −i ˆ 1 d ρˆ = [HS , ρ] ˆ + 2 γ˜ (ω)(([Sˆω , ρˆ Sˆω† ] + [Sˆω ρ, ˆ Sˆω† ]) dt ¯ 2¯ {ω≥0} ˆ Sˆω ])), + e−¯βω ([Sˆω† , ρˆ Sˆω ] + [Sˆω† ρ,
(2.12)
where HˆS =
k |kk|
(2.13)
is the renormalized physical system Hamiltonian. The system operators Sˆω are the Fourier transforms of the time evolved operators: ˆ ˆ −i Hˆ t/¯ = eiωt Sˆω , Sˆ−ω = Sˆω† , ei H t/¯ Se {ω}
{ω = ( k − l )/¯} − Bohr frequencies, (2.14) where Sˆ is the system operator which couples it to the bath. The temperature dependent friction function γ˜ (ω) is defined via the fluctuation dissipation relation for the autocorrelation
function of the bath induced force acting on the system: ∞ eitω Fˆ (t)Fˆ β dt, γ˜ (−ω) = e−¯βω γ˜ (ω), γ˜ (ω) = −∞
(2.15) ˆ ˆ where Fˆ (t) = ei HB t/¯ Fˆ e−i HB t/¯ . The thermal equilibrium Gibbs state ρˆβ = Z −1 exp{−β HˆS } is readily seen by inspection, to be a stationary solution of the master equation (2.12). C. Master equation for surface scattering
Consider then the Hamiltonian for the vertical motion of the particle as given in Eq. (2.1). The physical system Hamiltonian is identified as pˆ 2 HˆS = + V (ˆz), 2M and the bath induced force operator is seen to be Fˆ =
N √ M cj xˆj .
(2.16)
(2.17)
j =1
Averaging of the autocorrelation function over the thermal bath Hamiltonian and then carrying out the Fourier transform gives the thermal friction function as γ˜ (ω) =
¯M[J (ω) − J (−ω)] . 1 − exp(−¯βω)
(2.18)
For Ohmic friction we then identify γ˜ (ω) =
2M¯γ ω . 1 − exp(−¯βω)
(2.19)
The critical step is then to introduce a pair of adjoint ˆ Aˆ † defined as operators A, i 1 ˆ V (ˆz) ◦ p] Aˆ = [V (ˆz) + 2 Mω0
(2.20)
ˆ and Aˆ ◦ Bˆ ≡ (Aˆ Bˆ + Bˆ A)/2 is defined as a symmetric product of two self-adjoint operators. The structure of the operators Aˆ † , Aˆ (2.20) is uniquely (up to a frequency scale ω0 ) determined by the form of the particle Hamiltonian and its particular coupling to the bath. It does not depend on the detailed structure of the discrete part of the energy spectrum, which is not crucial for the scattering regime. The frequency scale ω0 describes a relevant energy scale E0 = ¯ω0 whose physical meaning is yet to be determined. We now propose the following Markovian master equation which describes scattering of a particle in the presence of a surface environment −i 1 d ρˆ ˆ ρˆ Aˆ † ] + [Aˆ ρ, ≡ L(ρ) ˆ = [HˆS , ρ]+ ˆ γ˜ (ω0 )(([A, ˆ Aˆ † ]) dt ¯ 2¯2 ˆ + [Aˆ † ρ, ˆ + e−¯βω0 ([Aˆ † , ρˆ A] ˆ A])). (2.21) This master equation describes both damping of the system by transferring energy to the bath as well as pumping of energy from the bath to the system. It is therefore well defined for all surface temperatures. Similar to Refs. 23–26 it is equivalent to a stochastic Ito-Schrödinger equation driven by two independent Brownian or Poisson processes as shown in Refs. 47 and 48 (a similar formalism was independently developed and
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-4
Li et al.
J. Chem. Phys. 140, 014104 (2014)
used for numerical simulations in Ref. 49 and many following papers). The semi-phenomenological master equation (2.21) possesses several properties which justify its applicability at least on time scales comparable or longer than 1/ω0 : (I) The solution of Eq. (2.21) is given by a completely positive, trace preserving semigroup. (II) For the case of a harmonic oscillator with V (ˆz) = Mω02 zˆ 2 /2 Eq. (2.21) coincides with the standard master equation for a harmonic oscillator linearly coupled to a heat bath at arbitrary temperature and satisfies all previous conditions (i)–(iv) enumerated above for the LGKS equation. (III) The dissipative part of the corrected (completely positive) Caldeira-Leggett equation can be obtained using (2.21) with a harmonic potential V (ˆz) and ω0 = kB T/¯ (the inverse of the so-called “thermal memory time”). (IV) The equation of motion obtained in the classical limit at T = 0 K, coincides with the Langevin equation given in Eq. (2.8) in the same limit (see below). Equation (2.21) is the central object of this paper. We will refer to it as the Space Dependent Friction Master Equation (SDFME), as it generalizes the standard master equation to the case in which the friction depends on the coordinate of the system. We note that at this point, it still depends on an arbitrary parameter, the frequency ω0 . This dependence is reminiscent of a similar free parameter – the width parameter – which is arbitrary when using frozen Gaussian semiclassical initial value representations for the quantum propagator.50 D. The energy balance for the Master equation
The time derivative of the kinetic energy observable Eˆ k = pˆ 2 /2M in the Heisenberg picture may be derived from the SDFME (Eq. (2.21)): d Eˆ k 1 ˆ + [Aˆ † , Eˆ k ]Aˆ = L∗ (Eˆ k ) = 2 γ˜ (ω0 )(Aˆ † [Eˆ k , A] dt 2¯ ˆ Eˆ k , Aˆ † ] + [A, ˆ Eˆ k ]Aˆ † )) (2.22) + e−¯βω0 (A[ and can be represented in terms of a temperature-independent damping power operator Pˆd Mγ ω0 ˆ † ˆ ˆ ˆ Pˆd = − (A [Ek , A] + [Aˆ † , Eˆ k ]A) ¯ and pumping power operator Pp
(2.23)
Mγ ω0 ˆ ˆ ˆ † ˆ Eˆ k ]Aˆ † ) − Pˆd (A[Ek , A ] + [A, Pˆp = ¯ Mγ ω0 ˆ ˆ ˆ † ˆ [A, [Ek , A ]] + [Aˆ † , [Eˆ k , A]]. = (2.24) ¯ The change of the average kinetic energy Eˆ k of the system is then obtained from the master equation for the kinetic energy of the particle as d dEˆ k = Tr(ρ(t) ˆ Eˆ k ) = Tr(ρ(t) ˆ L∗ (Eˆ k )) dt dt 1 Tr(ρ(t) ˆ Pˆp ). (2.25) = −Tr(ρ(t) ˆ Pˆd ) + ¯βω e 0 −1
E. Classical limit of the averaged dynamics
It is instructive to consider the averaged dynamics of the SDFME in the classical limit, when ¯ → 0. For this purpose one has first to compute the second time derivative of the averaged position Z = Tr(ρˆ zˆ ). For a mean value W = Tr(ρˆ w) ˆ of any observable wˆ we have the formula for the nth time derivative dn W = Tr(ρ[L ˆ ∗ ]n (w)). ˆ dt n
(2.26)
ˆ as the Heisenberg picture conjugate of the Denoting L∗ (ρ) Schrödinger picture generator L(ρ) ˆ defined in the SDFME (Eq. (2.21)) we have for the coordinate operator (for the sake of brevity we left out the “hat” notation) that L∗ (z) =
p 1 − γ˜ (ω0 )[1 − e−¯βω0 ] V V M 4¯ω0 M + γ˜ (ω0 )[1 + e−¯βω0 ]
1 V V , (2.27) 4(2ω0 M)2
and for the momentum: L∗ (p) = −V − γ˜ (ω0 )[1−e−¯βω0 ] − γ˜ (ω0 )[1+e−¯βω0 ]
1 [(V 2 −V V ] ◦ p 4¯ω0 M
1 [(V 2 −V V ] ◦ p. 4(2ω0 M)2 (2.28)
Similarly one can compute [L∗ ]2 (z). Keeping the lowest order terms in ¯ and using the Ehrenfest theorem for localized states, i.e., replacing the averages of functions of z by functions of the average Z = Tr(ρz) one obtains the Newton equation M
d d2 Z = −V (Z) − (Z) Z, 2 dt dt 1
(Z) = γ˜ (ω0 )[1 − e−¯βω0 ] V 2 . 2¯ω0
(2.29)
Notice that (a) Due to Eq. (2.18) (Z) does not depend on the temperature. (b) For Ohmic friction, from Eq. (2.19) we find that (Z) = Mγ V 2 , it is independent of ¯. (c) As expected, we have derived the averaged weak damping version of the classical Langevin equation, given in Eq. (2.8) or equivalently the T = 0 K limit for the classical dynamics. It is also of interest to consider explicitly the classical limit for the change in the kinetic energy of the particle. Keeping lowest order terms in ¯ for the damping and pumping power operator (Eqs. (2.23) and (2.24)) one has that γ (pˆ ◦ Vˆ )2 Pˆd = M
(2.30)
Pˆp = ¯γ ω0 (Vˆ )2 .
(2.31)
and
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-5
Li et al.
J. Chem. Phys. 140, 014104 (2014)
Noting that (e¯βω0 − 1)−1 (¯βω0 )−1 we find that the average change in the kinetic energy reduces in this limit to
18 17.8 17.6
ΔΕ (meV)
γ dE − Tr(ρ(t)( ˆ pˆ ◦ Vˆ )2 ) + γ kB T Tr(ρ(t)( ˆ Vˆ )2 ) dt M (2.32) and this has the same form as the classical result (Eq. (2.11)).
18.2
17.4 17.2 17 16.8
III. NUMERICAL RESULTS
16.6
As mentioned in the Introduction, the model we used is based on our previous fitting for the classical scattering of Ar on a LiF surface. The model parameters may thus be considered to be realistic and typical of atom surface scattering. For the vertical system motion we use a Morse potential (3.1)
with a well depth of V0 = 88 meV and a stiffness parameter α = 0.75 Å−1 . The mass of Ar is 39.948 a.m.u. Using Ohmic friction one defines a reduced friction parameter as √ 3/2 γ¯ = 2γ 2Mα 3 V0 (3.2) in all computations we then chose the value γ¯ = 3.38 × 10−3 or equivalently γ = 0.3851647229 a.u. The solution of the reduced quantum equation of motion was obtained by using the Fourier basis, Discrete Variable Representation (DVR) method, or sin-DVR method.51 The equally spaced spatial grid was chosen such that the minimum value of the vertical coordinate was −4.45 a.u.. The eigenfunctions associated with such a uniform grid are particle in a box eigenfunctions making it relatively easy to obtain the necessary matrix representation of the potential. The initial wavefunction was chosen to be a Gaussian such that Mω Mω 1/4 (z − z0 )2 exp − ψ(z, t = 0) = π¯ 2¯ i + pz0 (z − z0 ) , (3.3) ¯ with z0 = 32.0 a.u. and the width of the Gaussian determined by the frequency ω = 10−6 a.u. With this choice the spatial )−1/2 = 2.62 a.u., similar width of the initial density is ( 2Mω ¯ to the spatial width of the Morse potential measured as α −1 = 2.52 a.u. The initial momentum pz0 is determined by the initial kinetic energy of the particle. To implement the SDFME one needs to determine the unknown frequency parameter ω0 . We know that the exact quantum dynamics is independent of the free parameter. We therefore choose it as the value at which the derivative of the energy loss vanishes with respect to ω0 . In Fig. 1 we plot the numerically determined energy loss as found from the SDFME as a function of the frequency parameter ω0 for two different surface temperatures. We find a bell shaped dependence and a maximum which moves to higher frequencies as the temperature is increased. The extremal value (ω0∗ ) is then plotted in Fig. 2 as a function of the square root of the incident energy for two different
16.2
0
0.0005
0.001
0.0015
0.002
0.0025
ω 0 (a.u.)
FIG. 1. The dependence of the energy loss on the frequency parameter ω0 appearing in the SDFME. Results are shown for the surface temperatures T = 0 K (solid line) and T = 100 K (dashed line) at an incident energy of 100 meV. Note the single maximum in the plot.
surface temperatures. In both cases,√ a linear increasing dependence of the optimal frequency on Ei is found. This is reasonable, as ω0 sets the typical time scale of the collision. As the energy is increased the collision time decreases roughly inversely with the velocity so that one expects the frequency to increase linearly with the velocity. It is of interest to compare the quantum energy loss and its dependence on the incident energy and surface temperature with the classical result. In Fig. 3 we present the dependence of the average energy loss on the surface temperature, as obtained from the SDFME, for four different initial energies (solid line in all four panels). The dashed line in all four panels shows the classical energy loss predicted from perturbation theory, as given in Eq. (2.11). Since the momentum distribution of the initial wavepacket is also a Gaussian, with )1/2 (see also Eq. (3.3)) we computed momentum width ( M¯ω 2 classically the energy loss as obtained with such a distribution of the initial momenta centered about the average initial
0.001
0.0009
0.0008
ω0∗(a.u.)
V (z) = V0 [1 − exp(−αz)]2 − V0 ,
16.4
0.0007
0.0006
0.0005
0.0004
6
8
10
12
14
16
18
20
22
24
√Ei(meV) ⎯⎯
FIG. 2. Energy dependence of the optimized frequency ω0∗ . The optimized frequency parameter of the SDFME is plotted as a function of the square root of the incident energy, in the range 50 ≤ Ei ≤ 500 meV at the surface temperatures T = 0 K (solid line) and T = 100 K (dashed line).
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-6
Li et al.
J. Chem. Phys. 140, 014104 (2014)
(a)
(b) 85
ΔΕ(meV)
ΔΕ(meV)
195
175
155
75
65 0
50
100
150
200
250
300
0
50
100
T(K)
150
200
250
300
250
300
T(K)
(c)
(d) 12
ΔΕ(meV)
ΔΕ(meV)
20
18
10.5
9
16 0
50
100
150
200
250
300
0
50
100
T(K)
150
200
T(K)
FIG. 3. The temperature dependence of the energy loss. Panels (a)–(d) are at the incident energies of 500, 300, 100, and 50 meV, respectively. The solid, dotted, and dashed lines show the quantum results as obtained from the SDFME, the classical Langevin equation, and the leading order classical perturbation theory results.
the Ar atom increases to 0.27 a.u. the differences are not yet very large. However, the quantum effect, perhaps as expected increases as the energy is lowered. This is shown in Fig. 4 where we plot the relative change: μ(T = 0) =
ECl − EQ EQ
(3.4)
found at T = 0 K as a function of the incident energy.
0.08 0.075 0.07 0.065 0.06 0.055
μ
momentum. The dotted line in all four panels shows the numerically exact classical energy loss (obtained via numerical solution of the classical Langevin equation, using a second order method as described in Ref. 52 and averaged over the initial momentum distribution). The noticeable differences between the analytical estimate of the classical energy loss (Eq. (2.11)) and the numerically exact classical energy loss indicate that the friction employed is not really very weak. The analytic result should be considered more as a qualitative guide rather than a quantitative estimate. We also note that when T = 0 K, the quantum and classical energy losses are rather close to each other. Semiclassically it has been shown that to leading order in the friction coefficient the zero temperature semiclassical and classical energy losses should be identical.54 However, as already discussed, the numerically exact classical energy loss as obtained from the classical Langevin equation and the analytical energy loss as found from Eq. (2.11) are quite different for the parameters chosen for our simulation. The proximity of the quantum and classical T = 0 K energy losses is thus not self evident. As the temperature increases, as expected, the energy loss decreases both in the quantum and the classical simulations. At high incident energy (Ei = 500 meV, panel (a) of Fig. 3) the differences between the classical and quantum energy losses are rather small. At this energy the deBroglie wavelength of the Ar atom is 0.12 a.u., much smaller than the spatial extent of the Morse potential measured as 1/α = 2.52 a.u. Even at an energy of 100 meV, for which the deBroglie wavelength of
0.05 0.045 0.04 0.035 0.03 0.025 50
100
150
200
250
300
350
400
450
500
Ε i(meV)
FIG. 4. The incident energy dependence of the magnitude of the quantum energy loss as compared to the numerically exact classical energy loss (dotted lines in Fig. 3) at T = 0 K. For further details, see the text.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-7
Li et al.
J. Chem. Phys. 140, 014104 (2014)
140
The Heidelberg MCTDH Package32 was employed to compute the T = 0 K quantum dynamics of the Hamiltonian system as defined in Eq. (2.1), using the Morse potential for the system (Eq. (3.1)). The MCTDH wave function at time t = 0 is a Hartree product of one-dimensional functions:
130
(0) = ϕ(z)
N
ϕ(xj ).
(4.1)
j =1
ϕ(z) = N e
e
i ¯ pz0 (z−z0 )
.
ϕ(xj ) = Nj e− 2¯ ωj xj . 2
cj = ωj
100
2γ ωc . π (N + 1)
80
0
1000
2000 3000 Time [fs]
4000
5000
FIG. 5. Time dependence of the average kinetic energy for an Ar atom coupled to a single (ground state) bath oscillator. The results are presented for ωc = 0.0005 a.u. the value which gives the maximal energy loss.
(4.3)
The parameters of the system Hamiltonian as well as the friction coefficient were identical to those chosen for the simulation of the master equation (see Sec. III). The bath modes were discretized (with N + 1 modes) to mimic an Ohmic spectral density with a cutoff frequency ωc , using the discretization algorithm of Ref. 53 (assuming that the frequency of the N + 1-th mode is infinite):
j (4.4) ωj = −ωc log 1 − N +1 and
110
(4.2)
The initial wavefunction for the jth bath mode (using mass weighted coordinates) is its ground state wavefunction: 1
120
90
The system wavefunction is a Gaussian centered about z0 with incident average momentum pz0 : 1 − 2¯ Mω(z−z0 )2
[meV]
IV. COMPARISON WITH NUMERICALLY EXACT QUANTUM RESULTS
(4.5)
We employed the constant mean-field (CMF) integration scheme with a variable step size and an error tolerance of 10−8 . Within the CMF scheme the Bulirsch-Stoer (BS) integrator was employed for propagating the single particle functions (spf) and the short iterative Lanczos (SIL) integrator was employed for propagating the MCTDH coefficients. The error tolerance for the BS and SIL integrators was taken to be 10−8 . Provided that the integration scheme is accurate, the convergence of the MCTDH calculation to the numerically exact solution depends on the size of the DVR basis (grid points) and the number of single particle functions associated with each degree of freedom. In particular, the vertical (z) degree of freedom was described by a Fourier basis DVR and the bath (xj ) degrees of freedom were described by a Harmonic Oscillator (HO) DVR. The Fourier DVR was taken as 2025 grid points with the end points at z = −2.7 and z = 45 a.u. The HO DVR was taken as 90 grid points for each bath degree of freedom. The equilibrium position of each HO DVR was 0.0 and its frequency was determined according to Eq. (4.4). This choice is probably an overkill, intended to assure that the numerics is accurate. For the master equation the friction is purely Ohmic (the cutoff frequency ωc → ∞). For the MCTDH computation a cutoff frequency was introduced due to the fact that we are using a finite number of bath modes. To mimic the SDFME for
a fixed number of bath modes, we varied the cutoff frequency using the value which gave a maximal energy loss, this procedure is further discussed below. Here we note that Ohmic friction is expected to give a larger energy loss than a memory friction, justifying using the maximal average energy loss obtained by varying the cutoff frequency. The quantum dynamics was computed for Hamiltonians with N = 1, 2, 4, 8, and 12 modes. The time interval for the evolution of the wave function was 4998 fs, assuring that after the scattering event, the particle was essentially in the asymptotic region of the surface. Data from the calculation were registered to the output file every 17 fs. For the case of a single bath mode, the vertical degree of freedom was assigned with 20 single particle functions and the single mode was assigned with 20 single particle functions. A cutoff frequency of 0.0005 a.u. gives a maximum energy loss value of 12.1 meV. The time dependence of the average kinetic energy is plotted in Fig. 5. The energy loss values as a function of the cutoff frequencies are summarized in Table I. Note the bell shaped dependence of the energy loss on the cutoff frequency. For the simulation with two bath modes, the vertical degree of freedom was assigned with 32 single particle functions and the two bath modes were treated as a combined mode assigned with 32 single particle functions. The time dependence of the average kinetic energy is qualitatively similar to that found for the single bath mode case. A cutoff frequency of TABLE I. Cutoff frequency dependence of the average energy loss for a single bath mode. ωc (a.u.) 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009
E (meV) 5.9 9.9 12.1 11.9 10.1 7.8 5.5
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
Li et al.
J. Chem. Phys. 140, 014104 (2014)
TABLE II. Cutoff frequency dependence of the average energy loss for two bath modes.
0.0005 0.0007 0.0009 0.0011 0.0013
E (meV) 9.7 12.5 14.0 12.7 9.8
0.0009 a.u. gives a maximum energy loss value of 14.0 meV, as summarized in Table II. For 4 modes, the vertical degree of freedom was assigned with 32 single particle functions and modes {1, 2} and {3, 4} were treated as a combined mode assigned with 32 and 9 single particle functions, respectively. The cutoff frequency of 0.0016 a.u. gives a maximum energy loss value of 15.7 meV. The energy loss values as a function of the cutoff frequencies are summarized in Table III. For 8 modes, the vertical degree of freedom was assigned as before with 32 single particle functions and modes {1, 2}, {3, 4}, {5, 6}, and {7, 8} were treated as a combined mode assigned with 32, 8, 6, and 4 single particle functions, respectively. The cutoff frequency of 0.0028 a.u. gave a maximum energy loss value of 17.1 meV as summarized in Table IV. Finally the simulation with 12 bath modes was carried out with the vertical degree of freedom assigned with 26 single particle functions and modes {1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, and {11, 12} were treated as a combined mode assigned with 25, 7, 4, 4, 3, and 3 single particle functions, respectively. The maximal energy loss 17.7 meV was found with a cutoff frequency of 0.0040 a.u., see Table V. As the number of bath modes increases we note that the dependence of the average energy loss on the cutoff frequency decreases. The ratio of the curvatures for Tables I–V is 84:25:7.8:2.1:1. We also find that the maximal energy loss converges as a function of the number of bath modes. This convergence is shown in Fig. 6. Having established that the MCTDH computations are essentially converged when using 12 bath modes, we compared the results with those obtained at the same incident energy of 100 meV and T = 0 K from the SDFME computations presented in Sec. III. In Fig. 7 we compare the results for the time dependence of the average kinetic energy of the particle obtained with the two methods and find almost perfect agreement. TABLE III. Cutoff frequency dependence of the average energy loss for four bath modes.
E (meV)
ωc (a.u.) 0.0020 0.0022 0.0024 0.0026 0.0028 0.0030 0.0032
16.1 16.5 16.8 17.1 17.1 17.0 16.6
TABLE V. Cutoff frequency dependence of the average energy loss for 12 bath modes. E (meV)
ωc (a.u.) 0.0034 0.0036 0.0038 0.0040 0.0042 0.0044 0.0046
17.4 17.6 17.7 17.7 17.6 17.5 17.2
18 17
[meV]
ωc (a.u.)
TABLE IV. Cutoff frequency dependence of the average energy loss for eight bath modes.
16 15 14 13 12
1
2
4
8
12
N FIG. 6. The maximal average energy loss is plotted as a function of the number of bath modes (N) used in the simulation.
140
MCTDH Master Eq
130
[meV]
014104-8
120 110 100 90
ωc (a.u.) 0.0010 0.0012 0.0014 0.0016 0.0018 0.0020
E (meV) 13.6 14.7 15.5 15.7 15.1 13.8
80
0
1000
2000 3000 Time [fs]
4000
FIG. 7. Time dependence of the average kinetic energy as obtained from the SDFME (crosses) and MCTDH (plus signs) computations at an incident energy of 100 meV. The good agreement supports the validity of the SDFME approximation.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
Li et al.
014104-9
J. Chem. Phys. 140, 014104 (2014)
MCTDH Master Eq
30
[a.u.]
23
16
9
2 0
1000
2000 3000 Time [fs]
4000
FIG. 8. Time dependence of the average position of the incident Ar atom as obtained from the SDFME (crosses) and MCTDH (plus signs) computations for an incident energy of 100 meV. The good agreement supports the validity of the SDFME approximation.
The same good agreement is also found for the time dependence of the average position of the incident particle, as shown in Fig. 8. Finally in Fig. 9 we compare the time dependence of the diagonal matrix element of the reduced density operator as obtained from the MCTDH and SDFME computations at different times during the evolution of the scattering of the particle on the surface. Here too, the agreement is quantitative. V. DISCUSSION
A generalized Lindblad like master equation was written down for the scattering of an atom from a surface. The space dependent friction master equation has the necessary physical properties, especially including the nonlinear dependence of the interaction of the particle with the surface. The SDFME adheres to some general physical principles, such as being completely positive, trace preserving, and reducing to the standard master equation in the limit of a harmonic oscillator interaction potential.
[a.u.]
0.4
(a)
Time=0.0 [a.u.] Time=90000.0 [a.u.] Time=120000.0 [a.u.] Time=150000.0 [a.u.] Time=195000.0 [a.u.]
0.2
ACKNOWLEDGMENTS
0.1 0 0.4
[a.u.]
SDFME
0.3
However, the SDFME introduces an arbitrary time scale parameter, the frequency ω0 . Since the “true” quantum dynamics of the particle is independent of this parameter, it is chosen by demanding a stationarity of the solution for the relevant observable with respect to change of the frequency. We also found empirically, that the time scale determined by the frequency is linearly dependent on the velocity of the incident particle, thus reflecting the interaction time of the particle with the surface. One may expect this interaction time to be rather insensitive to the observable used for determining the optimal ω0 so that different observables will lead to similar frequencies so that in practice the SDFME is well defined. At vanishing temperature the average dynamics of the SDFME in the limit of ¯ → 0 reduces to the classical Langevin dynamics. At higher temperature, this remains correct only in the limit of weak damping. In this sense the SDFME differs from the Caldeira-Leggett master equation formalism. The Wigner transform of the SDFME in the limit ¯ → 0 is not identical to the corresponding Fokker-Planck equation. Numerical computations show that the differences are not big, however, it is important to note that they are not identical. The main reason for creating the SDFME is to study the quantum dynamics of the scattering. These will usually take place when the surface temperature is low. It was therefore important to show that indeed in the low temperature limit, the SDFME dynamics is “correct.” For this purpose, we also carried out numerically exact MCTDH computations for the Hamiltonian underlying the SDFME and found that at T = 0 K the agreement between the two dynamics is very good. In this paper, we have limited the discussion to scattering in the vertical direction only. The quantum effects measured were therefore not very large. Typically, the quantum energy loss at low temperature was lower than the classical, this may be related to the diffuseness of the quantum particle, leading to a weaker interaction with the surface. The more prominent quantum effects such as diffraction appear only if one adds the horizontal corrugation potential, that is one considers the SDFME for two or three degrees of freedom. Even for three degrees of freedom the SDFME may be solved numerically, so that this should serve as a good model for studying the quantum dynamics of scattering from a corrugated thermal surface.
0.3
0
10
(b)
20 MCTDH
30
40 z [a.u.]
Time=0.0 [a.u.] Time=90000.0 [a.u.] Time=120000.0 [a.u.] Time=150000.0 [a.u.] Time=195000.0 [a.u.]
0.2 0.1 0
0
10
20
30
40 z [a.u.]
This work was initiated while Robert Alicki was the Weston Visiting Professor at the Weizmann Institute of Science. It has been supported by grants from the National Natural Science Foundation of China, the 973 program of the Ministry of Science and Technology of China, the Israel Science Foundation, and the German-Israeli Foundation for Basic Research. 1 G.
Benedek and J. P. Toennies, Surf. Sci. 299, 587 (1994). Gross, Surf. Sci. Rep. 32, 291 (1998). 3 D. Farias and K. H. Rieder, Rep. Prog. Phys. 61, 1575 (1998). 4 B. Gumhalter, Phys. Rep. 351, 1 (2001). 5 A. S. Sanz and S. Miret-Artés, Phys. Rep. 451, 37 (2007). 2 A.
FIG. 9. Time dependence of the diagonal matrix element of the reduced density operator obtained by the master equation (panel (a)) and the MCTDH (panel (b)) at different propagation times at an incident energy of 100 meV.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-10
Li et al.
6 J. R. Manson, in Handbook of Surface Science, Energy Transfer to Phonons
in Atom and Molecule Collisions with Surfaces, edited by S. Holloway and N. V. Richardson (Elsevier, 2008), Vol. 3, p. 54. 7 A. T. Yinnon and R. Kosloff, Chem. Phys. Lett. 102, 216 (1983). 8 R. Kosloff, J. Phys. Chem. 92, 2087 (1988). 9 J. M. Moix, E. Pollak and W. Allison, J. Chem. Phys. 134, 024319 (2011). 10 H.-D. Meyer, Wiley Interdiscip. Rev.: Computat. Mol. Sci. 2, 351 (2012). 11 J. T. Stockburger, Chem. Phys. 296, 159 (2004). 12 Y. Zhou and J. Shao, J. Chem. Phys. 128, 034106 (2008). 13 H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003). 14 J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 170407 (2002). 15 J. Shao, J. Chem. Phys. 120, 5053 (2004). 16 S. Miret Artés and E. Pollak, Surf. Sci. Rep. 67, 161 (2012). 17 J. R. Manson, Comput. Phys. Commun. 80, 145 (1994). 18 T. Andersson, F. Althoff, P. Linde, S. Andersson, and K. Burke, Phys. Rev. B 65, 045409 (2002), and references therein. 19 J.M. Moix and E. Pollak, J. Chem. Phys. 134, 011103 (2011). 20 S. Miret-Artës, S. Daon and E. Pollak, J. Chem. Phys. 136, 204707 (2012). 21 K. Burke and W. Kohn, Phys. Rev. B 43, 2477 (1991). 22 E. Pollak, J. Phys. Chem. A 115, 7189 (2011), and references therein. 23 R. Alicki, Z. Phys. A 307, 279 (1982). 24 S. Gao, Phys. Rev. B 60, 15609 (1999). 25 S. Gao, J. Strømquist, and B. I. Lundquist, Phys. Rev. Lett. 86, 1805 (2001). 26 M. Nest, Phys. Rev. A 65, 052117 (2002). 27 V. Hakim and V. Ambegaokar, Phys. Rev. A 32, 423 1985). 28 B. Vacchini and K. Hornberger, Phys. Rep. 478, 71 (2009). 29 G. Lindblad, Commun. Math. Phys. 40, 147 (1975). 30 V. Gorini, A. Kossakowski and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976). 31 R. Alicki, Phys. Rev. A 40, 4077 (1989). 32 G. A. Worth, M. H. Beck, A. Jäckle, and H.-D. Meyer, The MCTDH Package, Version 8.2, University of Heidelberg, Heidelberg, Germany, 2000;
J. Chem. Phys. 140, 014104 (2014) H.-D. Meyer, Version 8.3, 2002; Version 8.4, 2007, see http://mctdh. uni-hd.de 33 H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990). 34 U. Manthe, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 97, 3199 (1992). 35 A. Capellini and A. P. J. Jansen, J. Chem. Phys. 104, 3366 (1996). 36 U. Manthe, J. Chem. Phys. 105, 6989 (1996). 37 M. Ehara, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 105, 8865 (1996). 38 M. H. Beck and H.-D. Meyer, Z. Phys. D 42, 113 (1997). 39 M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, Phys. Rep. 324, 1 (2000). 40 J. Trin, M. Monnerville, B. Pouilly, and H.-D. Meyer, J. Chem. Phys. 118, 600 (2003). 41 M. Nest and H.-D. Meyer, J. Chem. Phys. 119, 24 (2003). 42 H.-D. Meyer and G. A. Worth, Theor. Chem. Acc. 109, 251 (2003). 43 M. Nest and R. Kosloff, J. Chem. Phys. 127, 134711 (2007). 44 Multidimensional Quantum Dynamics: MCTDH Theory and Applications, edited by H.-D. Meyer, F. Gatti, and G. A. Worth (Wiley-VCH, Weinheim, 2009). 45 E. B. Davies, Commun. Math. Phys. 39, 91 (1974). 46 R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, 2nd ed., Lecture Notes in Physics Vol. 717 (Springer-Verlag, Berlin, 2007). 47 R. Alicki and M. Fannes, Lett. Math. Phys. 11, 259 (1986). 48 R. Alicki and M. Fannes, Commun. Math. Phys. 108, 353 (1987). 49 N. Gisin and I. C. Percival, J. Phys. A 25, 5677 (1992). 50 W. H. Miller, J. Chem. Phys. 125, 132305 (2006). 51 D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992). 52 E. V. Eijnden and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006). 53 H. Wang, M. Thoss and W. H. Miller, J. Chem. Phys. 112, 47 (2000). 54 I. Rips and E. Pollak, Phys. Rev. A 41, 5366 (1990).
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
THE JOURNAL OF CHEMICAL PHYSICS 140, 014104 (2014)
Quantum Markovian master equation for scattering from surfaces Haifeng Li,1 Jiushu Shao,1 Asaf Azuri,2 Eli Pollak,2,a) and Robert Alicki2,b) 1 2
College of Chemistry, Beijing Normal University, Beijing 100875, China Chemical Physics Department, Weizmann Institute of Science, 76100 Rehovoth, Israel
(Received 3 October 2013; accepted 30 November 2013; published online 6 January 2014) We propose a semi-phenomenological Markovian Master equation for describing the quantum dynamics of atom-surface scattering. It embodies the Lindblad-like structure and can describe both damping and pumping of energy between the system and the bath. It preserves positivity and correctly accounts for the vanishing of the interaction of the particle with the surface when the particle is distant from the surface. As a numerical test, we apply it to a model of an Ar atom scattered from a LiF surface, allowing for interaction only in the vertical direction. At low temperatures, we find that the quantum mechanical average energy loss is smaller than the classical energy loss. The numerical results obtained from the space dependent friction master equation are compared with numerical simulations for a discretized bath, using the multi-configurational time dependent Hartree methodology. The agreement between the two simulations is quantitative. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4851075] I. INTRODUCTION
The quantum dynamics of atom surface scattering presents a challenge to theory.1–6 It is relatively easy to solve for the quantum dynamics in the presence of a frozen surface. In this case, one has at most to undertake a three dimensional wavepacket propagation,7, 8 and this is routinely feasible, using, for example, split operator techniques.9 If however, one includes the interaction of the scattered particle with the surface, the problem becomes much more difficult. Even if the interacting surface modes are considered to be harmonic and the interaction is restricted to be linear in the coordinates of the surface modes, the problem is no longer straightforward due to the addition of the thermally distributed continuum of surface modes. In contrast to present day tractable problems10 such as the spin-boson model11, 12 or a system coupled bilinearly to a harmonic bath which may be treated using a variety of methods,13–15 here necessarily, the coupling is only linear in the bath coordinates but it must be nonlinear in the system coordinates. When the particle is far away from the surface the interaction vanishes and it is only felt when the particle is sufficiently close to the surface. It is this nonlinearity which makes a quantum solution for the dynamics difficult. One may well ask—who cares? Typically, if the particle is heavy, then the dynamics may be considered to be classical,16 while if it is sufficiently light, such as He one may use perturbation theory techniques6, 17 due to the very weak interaction of the atom with any surface. In fact though this is not sufficient. Quantum diffraction effects have been measured for heavy atoms and molecules scattered from surfaces.18 In most cases, the resulting angular distributions especially at very low surface temperatures have not received an adequate theoretical analysis. Moreover, it has been shown a) Electronic mail: [email protected] b) Permanent address: Institute of Theoretical Physics and Astrophysics,
University of Gda´nsk, Gda´nsk, Poland.
0021-9606/2014/140(1)/014104/10/$30.00
recently that one may observe quantum effects for heavy particle scattering simply by suitable collimation of the incident beam19, 20 provided that the de Broglie wavelength in the direction perpendicular to the motion of the incident particle is longer than the lattice length of the surface. There are also some good theoretical reasons for attempting to probe deeper into the quantum dynamics of surface scattering. For example, we do not know very much about the difference between the classical and quantum mechanical average energy lost by the particle when colliding with a surface.21 The same is true for the sticking probability to the surface as a function of incident energy or surface temperature.22 The central theme of this paper is to provide a framework for the study of these questions, based on a master equation description of the scattering dynamics. The formulation of a master equation for scattering is not a new problem. An early attempt to develop a formalism for dissipative scattering in terms of Markovian master equations was presented in the context of heavy-ion collisions in Ref. 23. A different approach was developed by Gao and co-workers24, 25 who introduce a master equation which is appropriate in the high-temperature regime as may be inferred by comparison with a generalization of the Caldeira-Leggett quantum Fokker-Planck type of equations. The theory used by Nest26 is of special relevance to the formulation presented in this paper. Using the formalism of super-symmetric quantum mechanics, Nest employed generalized raising and lowering operators for the system Hamiltonian. These operators encode the discrete structure of the Hamiltonian spectrum however they do not include the continuum part, which is inherent to the scattering process. Finally, we mention the class of master equations for a free particle interacting with a homogeneous equilibrium environment.27, 28 These equations satisfy the elementary conditions needed to describe the dissipative dynamics of a continuous spectrum Hamiltonian (see also below).
140, 014104-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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-2
Li et al.
In this paper we formulate a Lindblad like master equation, which introduces generalized Lindblad operators which are space dependent. They vanish when the particle is sufficiently far removed from the surface. In principle, the standard Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) form29, 30 is applicable under the condition that the motion of the system is relatively slow in comparison with the fast dynamics of a bath. The consistent derivation of such equations is available for the cases of well-separated time-scales, i.e., for spatially confined systems with a discrete spectrum of the Hamiltonian and in the weak coupling regime where dissipation is relatively slow in comparison with the system’s Hamiltonian motion. If this is not the case, different Markovian/nonMarkovian regimes are expected which depend on the time and energy scales determined by the initial conditions.31 As already noted, the scattering problems involve initial states (wave packets) composed of energy eigenstates from the continuous part of the system’s Hamiltonian spectrum. Therefore, we are not able to present a rigorous derivation but rather arguments based on extrapolation of the known results for the discrete spectrum case to our model and on the phenomenological properties of the solutions of master equations in the semiclassical regime. Since quantum effects will usually appear at low surface temperatures, it is important to verify that the theory presented is indeed valid in this regime. For this purpose, we compare the zero surface temperature results obtained from the master equation with a full quantum dynamics simulation, using a discretized bath and the Multi-Configurational Time Dependent Hartree (MCTDH) Heidelberg package.10, 32–44 We find excellent agreement, validating the usefulness of the master equation approach. In Sec. II we present the generalized quantum master equation. A numerical example, using parameters which are relevant for the scattering of an Ar atom from a LiF surface is presented in Sec. III. We find, that the average energy lost by the particle to the surface is especially at low temperatures, somewhat lower than the classical. In Sec. IV we describe the adaptation of the MCTDH methodology to the model studied with the master equation and compare the numerical results obtained at T = 0 K. We end with a discussion of the present methodology, its strengths, and weaknesses as well as its further generalization for treating the full three dimensional scattering problem.
II. QUANTUM MARKOVIAN MASTER EQUATIONS FOR DISSIPATIVE SCATTERING A. One-dimensional surface scattering model
In the following we will consider the dynamics of a particle with mass M incident on a smooth surface, interacting with it in the vertical direction only. This model may be readily generalized to include all three degrees of freedom (the vertical z and horizontal x and y directions), however the essential concepts are more readily demonstrated by considering only the vertical motion. The model Hamiltonian (with vertical momentum p) for describing the interaction of the system
J. Chem. Phys. 140, 014104 (2014)
with a harmonic bath takes the form16 ⎡ 2 ⎤ √ N Mc 1 ⎣ 2 p2 j +V (z)+ pj + ωj2 xj − H= V (z) ⎦ . 2M 2 j =1 ωj2 (2.1) The system potential V (z) vanishes for large vertical distances, it has the qualitative form of a Morse potential. The bath Hamiltonian using mass weighted coordinates (xj ) and momenta (pj ) is
1 2 p + ωj2 xj2 , 2 j =1 j N
HB =
(2.2)
while the system-bath interaction term is N √ HI = − MV (z) cj xj .
(2.3)
j =1
The continuum limit for the interaction with the surface is obtained by introducing the spectral density of bath modes: J (ω) =
N 2 π cj δ(ω − ωj ). 2 j =1 ωj
(2.4)
The friction function is then obtained through the Fourier transform of the spectral density 2 ∞ J (ω) γ (t) = . (2.5) dω cos(ωt) π 0 ω Since we want to approximate the dynamics with a Markovian master equation, we must consider the memoryless Ohmic spectral density J (ω) = γ ω
(2.6)
γ (t) = 2γ δ(t)
(2.7)
or conversely
identifying γ as the friction coefficient and δ(x) is the Dirac “delta” function. In the classical limit, with the Ohmic spectral density, the dynamics of the Hamiltonian given in Eq. (2.1) reduces to a space dependent Langevin equation of motion: √ (2.8) M z¨ + V (z) + Mγ [V (z)]2 z˙ = MV (z)F (t), where the (discretized) Gaussian random force is identified as
N pj F (t) = cj xj cos(ωj t) + sin(ωj t) . (2.9) ωj j =1 When averaged over the thermal bath distribution (exp(−βHB ) with β = 1/(kB T)) the random force vanishes, and its thermally averaged autocorrelation function is readily seen to satisfy the classical fluctuation dissipation relation: 1 F (t)F (t ) = γ (t − t ), (2.10) β with 1/β = kB T. Finally, we note that in the limit of weak friction, using classical perturbation theory one readily finds that the leading
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-3
Li et al.
J. Chem. Phys. 140, 014104 (2014)
order term for the average energy loss of the particle after colliding with the surface (E = Ei − Ef , where Ei and Ef are the initial and final energies of the scattered particle) is to first order in the friction ∞ 0 2 dV (zt ) γ ∞ 0 2 Eβ = Mγ dt − dt V zt , dt β −∞ −∞ (2.11) where zt0 denotes the solution of the unperturbed equation of motion M z¨ t0 + V (zt0 ) = 0 and ···β denotes the average with respect to the thermal state of the bath. The first term on the right hand side gives the average energy lost to the bath at T = 0 K, while the second term expresses the fact that as the temperature increases, the energy loss decreases. The bath can now also transfer energy back to the system.
B. LGKS equations for weakly coupled systems
Before writing down the reduced equation of motion for the scattering problem, we review briefly the structure of the Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) equation45, 46 in the weak coupling regime. The reduced equation of motion is constructed such that it preserves the following conditions: (i) The solution of the reduced equation of motion is given by a completely positive, trace preserving one-parameter semigroup. (ii) The diagonal elements of the density matrix in the energy basis evolve independently of the non-diagonal ones. (iii) The thermal canonical distribution is a stationary state. (iv) The detailed balance condition holds. The Schrödinger picture reduced dynamics of the system, in a weak coupling regime, is then given by the following Markovian master equation for the reduced density operator ρˆ −i ˆ 1 d ρˆ = [HS , ρ] ˆ + 2 γ˜ (ω)(([Sˆω , ρˆ Sˆω† ] + [Sˆω ρ, ˆ Sˆω† ]) dt ¯ 2¯ {ω≥0} ˆ Sˆω ])), + e−¯βω ([Sˆω† , ρˆ Sˆω ] + [Sˆω† ρ,
(2.12)
where HˆS =
k |kk|
(2.13)
is the renormalized physical system Hamiltonian. The system operators Sˆω are the Fourier transforms of the time evolved operators: ˆ ˆ −i Hˆ t/¯ = eiωt Sˆω , Sˆ−ω = Sˆω† , ei H t/¯ Se {ω}
{ω = ( k − l )/¯} − Bohr frequencies, (2.14) where Sˆ is the system operator which couples it to the bath. The temperature dependent friction function γ˜ (ω) is defined via the fluctuation dissipation relation for the autocorrelation
function of the bath induced force acting on the system: ∞ eitω Fˆ (t)Fˆ β dt, γ˜ (−ω) = e−¯βω γ˜ (ω), γ˜ (ω) = −∞
(2.15) ˆ ˆ where Fˆ (t) = ei HB t/¯ Fˆ e−i HB t/¯ . The thermal equilibrium Gibbs state ρˆβ = Z −1 exp{−β HˆS } is readily seen by inspection, to be a stationary solution of the master equation (2.12). C. Master equation for surface scattering
Consider then the Hamiltonian for the vertical motion of the particle as given in Eq. (2.1). The physical system Hamiltonian is identified as pˆ 2 HˆS = + V (ˆz), 2M and the bath induced force operator is seen to be Fˆ =
N √ M cj xˆj .
(2.16)
(2.17)
j =1
Averaging of the autocorrelation function over the thermal bath Hamiltonian and then carrying out the Fourier transform gives the thermal friction function as γ˜ (ω) =
¯M[J (ω) − J (−ω)] . 1 − exp(−¯βω)
(2.18)
For Ohmic friction we then identify γ˜ (ω) =
2M¯γ ω . 1 − exp(−¯βω)
(2.19)
The critical step is then to introduce a pair of adjoint ˆ Aˆ † defined as operators A, i 1 ˆ V (ˆz) ◦ p] Aˆ = [V (ˆz) + 2 Mω0
(2.20)
ˆ and Aˆ ◦ Bˆ ≡ (Aˆ Bˆ + Bˆ A)/2 is defined as a symmetric product of two self-adjoint operators. The structure of the operators Aˆ † , Aˆ (2.20) is uniquely (up to a frequency scale ω0 ) determined by the form of the particle Hamiltonian and its particular coupling to the bath. It does not depend on the detailed structure of the discrete part of the energy spectrum, which is not crucial for the scattering regime. The frequency scale ω0 describes a relevant energy scale E0 = ¯ω0 whose physical meaning is yet to be determined. We now propose the following Markovian master equation which describes scattering of a particle in the presence of a surface environment −i 1 d ρˆ ˆ ρˆ Aˆ † ] + [Aˆ ρ, ≡ L(ρ) ˆ = [HˆS , ρ]+ ˆ γ˜ (ω0 )(([A, ˆ Aˆ † ]) dt ¯ 2¯2 ˆ + [Aˆ † ρ, ˆ + e−¯βω0 ([Aˆ † , ρˆ A] ˆ A])). (2.21) This master equation describes both damping of the system by transferring energy to the bath as well as pumping of energy from the bath to the system. It is therefore well defined for all surface temperatures. Similar to Refs. 23–26 it is equivalent to a stochastic Ito-Schrödinger equation driven by two independent Brownian or Poisson processes as shown in Refs. 47 and 48 (a similar formalism was independently developed and
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-4
Li et al.
J. Chem. Phys. 140, 014104 (2014)
used for numerical simulations in Ref. 49 and many following papers). The semi-phenomenological master equation (2.21) possesses several properties which justify its applicability at least on time scales comparable or longer than 1/ω0 : (I) The solution of Eq. (2.21) is given by a completely positive, trace preserving semigroup. (II) For the case of a harmonic oscillator with V (ˆz) = Mω02 zˆ 2 /2 Eq. (2.21) coincides with the standard master equation for a harmonic oscillator linearly coupled to a heat bath at arbitrary temperature and satisfies all previous conditions (i)–(iv) enumerated above for the LGKS equation. (III) The dissipative part of the corrected (completely positive) Caldeira-Leggett equation can be obtained using (2.21) with a harmonic potential V (ˆz) and ω0 = kB T/¯ (the inverse of the so-called “thermal memory time”). (IV) The equation of motion obtained in the classical limit at T = 0 K, coincides with the Langevin equation given in Eq. (2.8) in the same limit (see below). Equation (2.21) is the central object of this paper. We will refer to it as the Space Dependent Friction Master Equation (SDFME), as it generalizes the standard master equation to the case in which the friction depends on the coordinate of the system. We note that at this point, it still depends on an arbitrary parameter, the frequency ω0 . This dependence is reminiscent of a similar free parameter – the width parameter – which is arbitrary when using frozen Gaussian semiclassical initial value representations for the quantum propagator.50 D. The energy balance for the Master equation
The time derivative of the kinetic energy observable Eˆ k = pˆ 2 /2M in the Heisenberg picture may be derived from the SDFME (Eq. (2.21)): d Eˆ k 1 ˆ + [Aˆ † , Eˆ k ]Aˆ = L∗ (Eˆ k ) = 2 γ˜ (ω0 )(Aˆ † [Eˆ k , A] dt 2¯ ˆ Eˆ k , Aˆ † ] + [A, ˆ Eˆ k ]Aˆ † )) (2.22) + e−¯βω0 (A[ and can be represented in terms of a temperature-independent damping power operator Pˆd Mγ ω0 ˆ † ˆ ˆ ˆ Pˆd = − (A [Ek , A] + [Aˆ † , Eˆ k ]A) ¯ and pumping power operator Pp
(2.23)
Mγ ω0 ˆ ˆ ˆ † ˆ Eˆ k ]Aˆ † ) − Pˆd (A[Ek , A ] + [A, Pˆp = ¯ Mγ ω0 ˆ ˆ ˆ † ˆ [A, [Ek , A ]] + [Aˆ † , [Eˆ k , A]]. = (2.24) ¯ The change of the average kinetic energy Eˆ k of the system is then obtained from the master equation for the kinetic energy of the particle as d dEˆ k = Tr(ρ(t) ˆ Eˆ k ) = Tr(ρ(t) ˆ L∗ (Eˆ k )) dt dt 1 Tr(ρ(t) ˆ Pˆp ). (2.25) = −Tr(ρ(t) ˆ Pˆd ) + ¯βω e 0 −1
E. Classical limit of the averaged dynamics
It is instructive to consider the averaged dynamics of the SDFME in the classical limit, when ¯ → 0. For this purpose one has first to compute the second time derivative of the averaged position Z = Tr(ρˆ zˆ ). For a mean value W = Tr(ρˆ w) ˆ of any observable wˆ we have the formula for the nth time derivative dn W = Tr(ρ[L ˆ ∗ ]n (w)). ˆ dt n
(2.26)
ˆ as the Heisenberg picture conjugate of the Denoting L∗ (ρ) Schrödinger picture generator L(ρ) ˆ defined in the SDFME (Eq. (2.21)) we have for the coordinate operator (for the sake of brevity we left out the “hat” notation) that L∗ (z) =
p 1 − γ˜ (ω0 )[1 − e−¯βω0 ] V V M 4¯ω0 M + γ˜ (ω0 )[1 + e−¯βω0 ]
1 V V , (2.27) 4(2ω0 M)2
and for the momentum: L∗ (p) = −V − γ˜ (ω0 )[1−e−¯βω0 ] − γ˜ (ω0 )[1+e−¯βω0 ]
1 [(V 2 −V V ] ◦ p 4¯ω0 M
1 [(V 2 −V V ] ◦ p. 4(2ω0 M)2 (2.28)
Similarly one can compute [L∗ ]2 (z). Keeping the lowest order terms in ¯ and using the Ehrenfest theorem for localized states, i.e., replacing the averages of functions of z by functions of the average Z = Tr(ρz) one obtains the Newton equation M
d d2 Z = −V (Z) − (Z) Z, 2 dt dt 1
(Z) = γ˜ (ω0 )[1 − e−¯βω0 ] V 2 . 2¯ω0
(2.29)
Notice that (a) Due to Eq. (2.18) (Z) does not depend on the temperature. (b) For Ohmic friction, from Eq. (2.19) we find that (Z) = Mγ V 2 , it is independent of ¯. (c) As expected, we have derived the averaged weak damping version of the classical Langevin equation, given in Eq. (2.8) or equivalently the T = 0 K limit for the classical dynamics. It is also of interest to consider explicitly the classical limit for the change in the kinetic energy of the particle. Keeping lowest order terms in ¯ for the damping and pumping power operator (Eqs. (2.23) and (2.24)) one has that γ (pˆ ◦ Vˆ )2 Pˆd = M
(2.30)
Pˆp = ¯γ ω0 (Vˆ )2 .
(2.31)
and
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-5
Li et al.
J. Chem. Phys. 140, 014104 (2014)
Noting that (e¯βω0 − 1)−1 (¯βω0 )−1 we find that the average change in the kinetic energy reduces in this limit to
18 17.8 17.6
ΔΕ (meV)
γ dE − Tr(ρ(t)( ˆ pˆ ◦ Vˆ )2 ) + γ kB T Tr(ρ(t)( ˆ Vˆ )2 ) dt M (2.32) and this has the same form as the classical result (Eq. (2.11)).
18.2
17.4 17.2 17 16.8
III. NUMERICAL RESULTS
16.6
As mentioned in the Introduction, the model we used is based on our previous fitting for the classical scattering of Ar on a LiF surface. The model parameters may thus be considered to be realistic and typical of atom surface scattering. For the vertical system motion we use a Morse potential (3.1)
with a well depth of V0 = 88 meV and a stiffness parameter α = 0.75 Å−1 . The mass of Ar is 39.948 a.m.u. Using Ohmic friction one defines a reduced friction parameter as √ 3/2 γ¯ = 2γ 2Mα 3 V0 (3.2) in all computations we then chose the value γ¯ = 3.38 × 10−3 or equivalently γ = 0.3851647229 a.u. The solution of the reduced quantum equation of motion was obtained by using the Fourier basis, Discrete Variable Representation (DVR) method, or sin-DVR method.51 The equally spaced spatial grid was chosen such that the minimum value of the vertical coordinate was −4.45 a.u.. The eigenfunctions associated with such a uniform grid are particle in a box eigenfunctions making it relatively easy to obtain the necessary matrix representation of the potential. The initial wavefunction was chosen to be a Gaussian such that Mω Mω 1/4 (z − z0 )2 exp − ψ(z, t = 0) = π¯ 2¯ i + pz0 (z − z0 ) , (3.3) ¯ with z0 = 32.0 a.u. and the width of the Gaussian determined by the frequency ω = 10−6 a.u. With this choice the spatial )−1/2 = 2.62 a.u., similar width of the initial density is ( 2Mω ¯ to the spatial width of the Morse potential measured as α −1 = 2.52 a.u. The initial momentum pz0 is determined by the initial kinetic energy of the particle. To implement the SDFME one needs to determine the unknown frequency parameter ω0 . We know that the exact quantum dynamics is independent of the free parameter. We therefore choose it as the value at which the derivative of the energy loss vanishes with respect to ω0 . In Fig. 1 we plot the numerically determined energy loss as found from the SDFME as a function of the frequency parameter ω0 for two different surface temperatures. We find a bell shaped dependence and a maximum which moves to higher frequencies as the temperature is increased. The extremal value (ω0∗ ) is then plotted in Fig. 2 as a function of the square root of the incident energy for two different
16.2
0
0.0005
0.001
0.0015
0.002
0.0025
ω 0 (a.u.)
FIG. 1. The dependence of the energy loss on the frequency parameter ω0 appearing in the SDFME. Results are shown for the surface temperatures T = 0 K (solid line) and T = 100 K (dashed line) at an incident energy of 100 meV. Note the single maximum in the plot.
surface temperatures. In both cases,√ a linear increasing dependence of the optimal frequency on Ei is found. This is reasonable, as ω0 sets the typical time scale of the collision. As the energy is increased the collision time decreases roughly inversely with the velocity so that one expects the frequency to increase linearly with the velocity. It is of interest to compare the quantum energy loss and its dependence on the incident energy and surface temperature with the classical result. In Fig. 3 we present the dependence of the average energy loss on the surface temperature, as obtained from the SDFME, for four different initial energies (solid line in all four panels). The dashed line in all four panels shows the classical energy loss predicted from perturbation theory, as given in Eq. (2.11). Since the momentum distribution of the initial wavepacket is also a Gaussian, with )1/2 (see also Eq. (3.3)) we computed momentum width ( M¯ω 2 classically the energy loss as obtained with such a distribution of the initial momenta centered about the average initial
0.001
0.0009
0.0008
ω0∗(a.u.)
V (z) = V0 [1 − exp(−αz)]2 − V0 ,
16.4
0.0007
0.0006
0.0005
0.0004
6
8
10
12
14
16
18
20
22
24
√Ei(meV) ⎯⎯
FIG. 2. Energy dependence of the optimized frequency ω0∗ . The optimized frequency parameter of the SDFME is plotted as a function of the square root of the incident energy, in the range 50 ≤ Ei ≤ 500 meV at the surface temperatures T = 0 K (solid line) and T = 100 K (dashed line).
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-6
Li et al.
J. Chem. Phys. 140, 014104 (2014)
(a)
(b) 85
ΔΕ(meV)
ΔΕ(meV)
195
175
155
75
65 0
50
100
150
200
250
300
0
50
100
T(K)
150
200
250
300
250
300
T(K)
(c)
(d) 12
ΔΕ(meV)
ΔΕ(meV)
20
18
10.5
9
16 0
50
100
150
200
250
300
0
50
100
T(K)
150
200
T(K)
FIG. 3. The temperature dependence of the energy loss. Panels (a)–(d) are at the incident energies of 500, 300, 100, and 50 meV, respectively. The solid, dotted, and dashed lines show the quantum results as obtained from the SDFME, the classical Langevin equation, and the leading order classical perturbation theory results.
the Ar atom increases to 0.27 a.u. the differences are not yet very large. However, the quantum effect, perhaps as expected increases as the energy is lowered. This is shown in Fig. 4 where we plot the relative change: μ(T = 0) =
ECl − EQ EQ
(3.4)
found at T = 0 K as a function of the incident energy.
0.08 0.075 0.07 0.065 0.06 0.055
μ
momentum. The dotted line in all four panels shows the numerically exact classical energy loss (obtained via numerical solution of the classical Langevin equation, using a second order method as described in Ref. 52 and averaged over the initial momentum distribution). The noticeable differences between the analytical estimate of the classical energy loss (Eq. (2.11)) and the numerically exact classical energy loss indicate that the friction employed is not really very weak. The analytic result should be considered more as a qualitative guide rather than a quantitative estimate. We also note that when T = 0 K, the quantum and classical energy losses are rather close to each other. Semiclassically it has been shown that to leading order in the friction coefficient the zero temperature semiclassical and classical energy losses should be identical.54 However, as already discussed, the numerically exact classical energy loss as obtained from the classical Langevin equation and the analytical energy loss as found from Eq. (2.11) are quite different for the parameters chosen for our simulation. The proximity of the quantum and classical T = 0 K energy losses is thus not self evident. As the temperature increases, as expected, the energy loss decreases both in the quantum and the classical simulations. At high incident energy (Ei = 500 meV, panel (a) of Fig. 3) the differences between the classical and quantum energy losses are rather small. At this energy the deBroglie wavelength of the Ar atom is 0.12 a.u., much smaller than the spatial extent of the Morse potential measured as 1/α = 2.52 a.u. Even at an energy of 100 meV, for which the deBroglie wavelength of
0.05 0.045 0.04 0.035 0.03 0.025 50
100
150
200
250
300
350
400
450
500
Ε i(meV)
FIG. 4. The incident energy dependence of the magnitude of the quantum energy loss as compared to the numerically exact classical energy loss (dotted lines in Fig. 3) at T = 0 K. For further details, see the text.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-7
Li et al.
J. Chem. Phys. 140, 014104 (2014)
140
The Heidelberg MCTDH Package32 was employed to compute the T = 0 K quantum dynamics of the Hamiltonian system as defined in Eq. (2.1), using the Morse potential for the system (Eq. (3.1)). The MCTDH wave function at time t = 0 is a Hartree product of one-dimensional functions:
130
(0) = ϕ(z)
N
ϕ(xj ).
(4.1)
j =1
ϕ(z) = N e
e
i ¯ pz0 (z−z0 )
.
ϕ(xj ) = Nj e− 2¯ ωj xj . 2
cj = ωj
100
2γ ωc . π (N + 1)
80
0
1000
2000 3000 Time [fs]
4000
5000
FIG. 5. Time dependence of the average kinetic energy for an Ar atom coupled to a single (ground state) bath oscillator. The results are presented for ωc = 0.0005 a.u. the value which gives the maximal energy loss.
(4.3)
The parameters of the system Hamiltonian as well as the friction coefficient were identical to those chosen for the simulation of the master equation (see Sec. III). The bath modes were discretized (with N + 1 modes) to mimic an Ohmic spectral density with a cutoff frequency ωc , using the discretization algorithm of Ref. 53 (assuming that the frequency of the N + 1-th mode is infinite):
j (4.4) ωj = −ωc log 1 − N +1 and
110
(4.2)
The initial wavefunction for the jth bath mode (using mass weighted coordinates) is its ground state wavefunction: 1
120
90
The system wavefunction is a Gaussian centered about z0 with incident average momentum pz0 : 1 − 2¯ Mω(z−z0 )2
[meV]
IV. COMPARISON WITH NUMERICALLY EXACT QUANTUM RESULTS
(4.5)
We employed the constant mean-field (CMF) integration scheme with a variable step size and an error tolerance of 10−8 . Within the CMF scheme the Bulirsch-Stoer (BS) integrator was employed for propagating the single particle functions (spf) and the short iterative Lanczos (SIL) integrator was employed for propagating the MCTDH coefficients. The error tolerance for the BS and SIL integrators was taken to be 10−8 . Provided that the integration scheme is accurate, the convergence of the MCTDH calculation to the numerically exact solution depends on the size of the DVR basis (grid points) and the number of single particle functions associated with each degree of freedom. In particular, the vertical (z) degree of freedom was described by a Fourier basis DVR and the bath (xj ) degrees of freedom were described by a Harmonic Oscillator (HO) DVR. The Fourier DVR was taken as 2025 grid points with the end points at z = −2.7 and z = 45 a.u. The HO DVR was taken as 90 grid points for each bath degree of freedom. The equilibrium position of each HO DVR was 0.0 and its frequency was determined according to Eq. (4.4). This choice is probably an overkill, intended to assure that the numerics is accurate. For the master equation the friction is purely Ohmic (the cutoff frequency ωc → ∞). For the MCTDH computation a cutoff frequency was introduced due to the fact that we are using a finite number of bath modes. To mimic the SDFME for
a fixed number of bath modes, we varied the cutoff frequency using the value which gave a maximal energy loss, this procedure is further discussed below. Here we note that Ohmic friction is expected to give a larger energy loss than a memory friction, justifying using the maximal average energy loss obtained by varying the cutoff frequency. The quantum dynamics was computed for Hamiltonians with N = 1, 2, 4, 8, and 12 modes. The time interval for the evolution of the wave function was 4998 fs, assuring that after the scattering event, the particle was essentially in the asymptotic region of the surface. Data from the calculation were registered to the output file every 17 fs. For the case of a single bath mode, the vertical degree of freedom was assigned with 20 single particle functions and the single mode was assigned with 20 single particle functions. A cutoff frequency of 0.0005 a.u. gives a maximum energy loss value of 12.1 meV. The time dependence of the average kinetic energy is plotted in Fig. 5. The energy loss values as a function of the cutoff frequencies are summarized in Table I. Note the bell shaped dependence of the energy loss on the cutoff frequency. For the simulation with two bath modes, the vertical degree of freedom was assigned with 32 single particle functions and the two bath modes were treated as a combined mode assigned with 32 single particle functions. The time dependence of the average kinetic energy is qualitatively similar to that found for the single bath mode case. A cutoff frequency of TABLE I. Cutoff frequency dependence of the average energy loss for a single bath mode. ωc (a.u.) 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009
E (meV) 5.9 9.9 12.1 11.9 10.1 7.8 5.5
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
Li et al.
J. Chem. Phys. 140, 014104 (2014)
TABLE II. Cutoff frequency dependence of the average energy loss for two bath modes.
0.0005 0.0007 0.0009 0.0011 0.0013
E (meV) 9.7 12.5 14.0 12.7 9.8
0.0009 a.u. gives a maximum energy loss value of 14.0 meV, as summarized in Table II. For 4 modes, the vertical degree of freedom was assigned with 32 single particle functions and modes {1, 2} and {3, 4} were treated as a combined mode assigned with 32 and 9 single particle functions, respectively. The cutoff frequency of 0.0016 a.u. gives a maximum energy loss value of 15.7 meV. The energy loss values as a function of the cutoff frequencies are summarized in Table III. For 8 modes, the vertical degree of freedom was assigned as before with 32 single particle functions and modes {1, 2}, {3, 4}, {5, 6}, and {7, 8} were treated as a combined mode assigned with 32, 8, 6, and 4 single particle functions, respectively. The cutoff frequency of 0.0028 a.u. gave a maximum energy loss value of 17.1 meV as summarized in Table IV. Finally the simulation with 12 bath modes was carried out with the vertical degree of freedom assigned with 26 single particle functions and modes {1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, and {11, 12} were treated as a combined mode assigned with 25, 7, 4, 4, 3, and 3 single particle functions, respectively. The maximal energy loss 17.7 meV was found with a cutoff frequency of 0.0040 a.u., see Table V. As the number of bath modes increases we note that the dependence of the average energy loss on the cutoff frequency decreases. The ratio of the curvatures for Tables I–V is 84:25:7.8:2.1:1. We also find that the maximal energy loss converges as a function of the number of bath modes. This convergence is shown in Fig. 6. Having established that the MCTDH computations are essentially converged when using 12 bath modes, we compared the results with those obtained at the same incident energy of 100 meV and T = 0 K from the SDFME computations presented in Sec. III. In Fig. 7 we compare the results for the time dependence of the average kinetic energy of the particle obtained with the two methods and find almost perfect agreement. TABLE III. Cutoff frequency dependence of the average energy loss for four bath modes.
E (meV)
ωc (a.u.) 0.0020 0.0022 0.0024 0.0026 0.0028 0.0030 0.0032
16.1 16.5 16.8 17.1 17.1 17.0 16.6
TABLE V. Cutoff frequency dependence of the average energy loss for 12 bath modes. E (meV)
ωc (a.u.) 0.0034 0.0036 0.0038 0.0040 0.0042 0.0044 0.0046
17.4 17.6 17.7 17.7 17.6 17.5 17.2
18 17
[meV]
ωc (a.u.)
TABLE IV. Cutoff frequency dependence of the average energy loss for eight bath modes.
16 15 14 13 12
1
2
4
8
12
N FIG. 6. The maximal average energy loss is plotted as a function of the number of bath modes (N) used in the simulation.
140
MCTDH Master Eq
130
[meV]
014104-8
120 110 100 90
ωc (a.u.) 0.0010 0.0012 0.0014 0.0016 0.0018 0.0020
E (meV) 13.6 14.7 15.5 15.7 15.1 13.8
80
0
1000
2000 3000 Time [fs]
4000
FIG. 7. Time dependence of the average kinetic energy as obtained from the SDFME (crosses) and MCTDH (plus signs) computations at an incident energy of 100 meV. The good agreement supports the validity of the SDFME approximation.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
Li et al.
014104-9
J. Chem. Phys. 140, 014104 (2014)
MCTDH Master Eq
30
[a.u.]
23
16
9
2 0
1000
2000 3000 Time [fs]
4000
FIG. 8. Time dependence of the average position of the incident Ar atom as obtained from the SDFME (crosses) and MCTDH (plus signs) computations for an incident energy of 100 meV. The good agreement supports the validity of the SDFME approximation.
The same good agreement is also found for the time dependence of the average position of the incident particle, as shown in Fig. 8. Finally in Fig. 9 we compare the time dependence of the diagonal matrix element of the reduced density operator as obtained from the MCTDH and SDFME computations at different times during the evolution of the scattering of the particle on the surface. Here too, the agreement is quantitative. V. DISCUSSION
A generalized Lindblad like master equation was written down for the scattering of an atom from a surface. The space dependent friction master equation has the necessary physical properties, especially including the nonlinear dependence of the interaction of the particle with the surface. The SDFME adheres to some general physical principles, such as being completely positive, trace preserving, and reducing to the standard master equation in the limit of a harmonic oscillator interaction potential.
[a.u.]
0.4
(a)
Time=0.0 [a.u.] Time=90000.0 [a.u.] Time=120000.0 [a.u.] Time=150000.0 [a.u.] Time=195000.0 [a.u.]
0.2
ACKNOWLEDGMENTS
0.1 0 0.4
[a.u.]
SDFME
0.3
However, the SDFME introduces an arbitrary time scale parameter, the frequency ω0 . Since the “true” quantum dynamics of the particle is independent of this parameter, it is chosen by demanding a stationarity of the solution for the relevant observable with respect to change of the frequency. We also found empirically, that the time scale determined by the frequency is linearly dependent on the velocity of the incident particle, thus reflecting the interaction time of the particle with the surface. One may expect this interaction time to be rather insensitive to the observable used for determining the optimal ω0 so that different observables will lead to similar frequencies so that in practice the SDFME is well defined. At vanishing temperature the average dynamics of the SDFME in the limit of ¯ → 0 reduces to the classical Langevin dynamics. At higher temperature, this remains correct only in the limit of weak damping. In this sense the SDFME differs from the Caldeira-Leggett master equation formalism. The Wigner transform of the SDFME in the limit ¯ → 0 is not identical to the corresponding Fokker-Planck equation. Numerical computations show that the differences are not big, however, it is important to note that they are not identical. The main reason for creating the SDFME is to study the quantum dynamics of the scattering. These will usually take place when the surface temperature is low. It was therefore important to show that indeed in the low temperature limit, the SDFME dynamics is “correct.” For this purpose, we also carried out numerically exact MCTDH computations for the Hamiltonian underlying the SDFME and found that at T = 0 K the agreement between the two dynamics is very good. In this paper, we have limited the discussion to scattering in the vertical direction only. The quantum effects measured were therefore not very large. Typically, the quantum energy loss at low temperature was lower than the classical, this may be related to the diffuseness of the quantum particle, leading to a weaker interaction with the surface. The more prominent quantum effects such as diffraction appear only if one adds the horizontal corrugation potential, that is one considers the SDFME for two or three degrees of freedom. Even for three degrees of freedom the SDFME may be solved numerically, so that this should serve as a good model for studying the quantum dynamics of scattering from a corrugated thermal surface.
0.3
0
10
(b)
20 MCTDH
30
40 z [a.u.]
Time=0.0 [a.u.] Time=90000.0 [a.u.] Time=120000.0 [a.u.] Time=150000.0 [a.u.] Time=195000.0 [a.u.]
0.2 0.1 0
0
10
20
30
40 z [a.u.]
This work was initiated while Robert Alicki was the Weston Visiting Professor at the Weizmann Institute of Science. It has been supported by grants from the National Natural Science Foundation of China, the 973 program of the Ministry of Science and Technology of China, the Israel Science Foundation, and the German-Israeli Foundation for Basic Research. 1 G.
Benedek and J. P. Toennies, Surf. Sci. 299, 587 (1994). Gross, Surf. Sci. Rep. 32, 291 (1998). 3 D. Farias and K. H. Rieder, Rep. Prog. Phys. 61, 1575 (1998). 4 B. Gumhalter, Phys. Rep. 351, 1 (2001). 5 A. S. Sanz and S. Miret-Artés, Phys. Rep. 451, 37 (2007). 2 A.
FIG. 9. Time dependence of the diagonal matrix element of the reduced density operator obtained by the master equation (panel (a)) and the MCTDH (panel (b)) at different propagation times at an incident energy of 100 meV.
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
014104-10
Li et al.
6 J. R. Manson, in Handbook of Surface Science, Energy Transfer to Phonons
in Atom and Molecule Collisions with Surfaces, edited by S. Holloway and N. V. Richardson (Elsevier, 2008), Vol. 3, p. 54. 7 A. T. Yinnon and R. Kosloff, Chem. Phys. Lett. 102, 216 (1983). 8 R. Kosloff, J. Phys. Chem. 92, 2087 (1988). 9 J. M. Moix, E. Pollak and W. Allison, J. Chem. Phys. 134, 024319 (2011). 10 H.-D. Meyer, Wiley Interdiscip. Rev.: Computat. Mol. Sci. 2, 351 (2012). 11 J. T. Stockburger, Chem. Phys. 296, 159 (2004). 12 Y. Zhou and J. Shao, J. Chem. Phys. 128, 034106 (2008). 13 H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003). 14 J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 170407 (2002). 15 J. Shao, J. Chem. Phys. 120, 5053 (2004). 16 S. Miret Artés and E. Pollak, Surf. Sci. Rep. 67, 161 (2012). 17 J. R. Manson, Comput. Phys. Commun. 80, 145 (1994). 18 T. Andersson, F. Althoff, P. Linde, S. Andersson, and K. Burke, Phys. Rev. B 65, 045409 (2002), and references therein. 19 J.M. Moix and E. Pollak, J. Chem. Phys. 134, 011103 (2011). 20 S. Miret-Artës, S. Daon and E. Pollak, J. Chem. Phys. 136, 204707 (2012). 21 K. Burke and W. Kohn, Phys. Rev. B 43, 2477 (1991). 22 E. Pollak, J. Phys. Chem. A 115, 7189 (2011), and references therein. 23 R. Alicki, Z. Phys. A 307, 279 (1982). 24 S. Gao, Phys. Rev. B 60, 15609 (1999). 25 S. Gao, J. Strømquist, and B. I. Lundquist, Phys. Rev. Lett. 86, 1805 (2001). 26 M. Nest, Phys. Rev. A 65, 052117 (2002). 27 V. Hakim and V. Ambegaokar, Phys. Rev. A 32, 423 1985). 28 B. Vacchini and K. Hornberger, Phys. Rep. 478, 71 (2009). 29 G. Lindblad, Commun. Math. Phys. 40, 147 (1975). 30 V. Gorini, A. Kossakowski and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976). 31 R. Alicki, Phys. Rev. A 40, 4077 (1989). 32 G. A. Worth, M. H. Beck, A. Jäckle, and H.-D. Meyer, The MCTDH Package, Version 8.2, University of Heidelberg, Heidelberg, Germany, 2000;
J. Chem. Phys. 140, 014104 (2014) H.-D. Meyer, Version 8.3, 2002; Version 8.4, 2007, see http://mctdh. uni-hd.de 33 H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990). 34 U. Manthe, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 97, 3199 (1992). 35 A. Capellini and A. P. J. Jansen, J. Chem. Phys. 104, 3366 (1996). 36 U. Manthe, J. Chem. Phys. 105, 6989 (1996). 37 M. Ehara, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 105, 8865 (1996). 38 M. H. Beck and H.-D. Meyer, Z. Phys. D 42, 113 (1997). 39 M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, Phys. Rep. 324, 1 (2000). 40 J. Trin, M. Monnerville, B. Pouilly, and H.-D. Meyer, J. Chem. Phys. 118, 600 (2003). 41 M. Nest and H.-D. Meyer, J. Chem. Phys. 119, 24 (2003). 42 H.-D. Meyer and G. A. Worth, Theor. Chem. Acc. 109, 251 (2003). 43 M. Nest and R. Kosloff, J. Chem. Phys. 127, 134711 (2007). 44 Multidimensional Quantum Dynamics: MCTDH Theory and Applications, edited by H.-D. Meyer, F. Gatti, and G. A. Worth (Wiley-VCH, Weinheim, 2009). 45 E. B. Davies, Commun. Math. Phys. 39, 91 (1974). 46 R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, 2nd ed., Lecture Notes in Physics Vol. 717 (Springer-Verlag, Berlin, 2007). 47 R. Alicki and M. Fannes, Lett. Math. Phys. 11, 259 (1986). 48 R. Alicki and M. Fannes, Commun. Math. Phys. 108, 353 (1987). 49 N. Gisin and I. C. Percival, J. Phys. A 25, 5677 (1992). 50 W. H. Miller, J. Chem. Phys. 125, 132305 (2006). 51 D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992). 52 E. V. Eijnden and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006). 53 H. Wang, M. Thoss and W. H. Miller, J. Chem. Phys. 112, 47 (2000). 54 I. Rips and E. Pollak, Phys. Rev. A 41, 5366 (1990).
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: 137.149.200.5 On: Mon, 01 Dec 2014 06:29:21
