and to validate the surface-hopping method include multitrajectory meth-

ods that preserve the coherence between trajectories. Although not often

applied, the mixed quantum-classical system could also be treated using

the Bohmian particle description of quantum mechanics, introduced in Sec-

tion 3.4. A Bohmian trajectory is a sample from the possible evolutions, and

the selection of a particular trajectory allows to handle the back reaction

114 Dynamics of mixed quantum/classical systems

for that particular trajectory (Section 5.3.4). In order to be able to handle

the full back reaction, an ensemble of Bohmian trajectories will be needed.

The mixed quantum-classical behavior is not limited to electrons and nu-

clei; we can just as well treat the quantum behavior of selected nuclei (no-

tably the proton) in a dynamical classical environment. In this case the

electrons are treated in the full Born“Oppenheimer approximation; they are

consistently in the electronic ground state and provide a potential ¬eld for

the nuclear interactions. But the quantum nucleus has the ability to tunnel

between energy wells, and move non-adiabatically with the involvement of

nuclear excited states. It should be noted that the quantum e¬ects of nu-

clei, being so much closer to classical behavior than electrons, can often be

treated by the use of e¬ective potentials in classical simulations, or even by

quantum corrections to classical behavior. Section 3.5 is devoted to these

approximate methods.

5.2 Quantum dynamics in a non-stationary potential

Assume that a quantum system of n particles r 1 , . . .r n interacts with a

time-dependent Hamiltonian

ˆ ˆ ˆ

H(t) = K + V (t). (5.3)

Then the time-dependent Schr¨dinger equation

o

‚ iˆ

Ψ(r, t) = ’ H(t)Ψ(r, t) (5.4)

‚t

can formally be solved as (see Chapter 14 for details of exponential opera-

tors)

t

i ˆ

Ψ(r, t) = exp ’ H(t ) dt Ψ(r, 0). (5.5)

0

This way of writing helps us to be concise, but does not help very much to

actually solve the time-dependent wave function.

In most applications a time-dependent solution in terms of expansion on

a basis set is more suitable than a direct consideration of the wave function

itself. If the time-dependence is weak, i.e., if the time-dependent interaction

can be considered as a perturbation, or if the time dependence arises from

the parametric dependence on slowly varying nuclear coordinates, the solu-

tion approaches a steady-state. Using expansion in a set of eigenfunctions

of the time-independent Schr¨dinger equation, all the hard work is then

o

done separately, and the time-dependence is expressed as the way the time-

independent solutions mix as a result of the time-dependent perturbation.

5.2 Quantum dynamics in a non-stationary potential 115

The basis functions in which the time-dependent wave function is ex-

panded can be either stationary or time-dependent, depending on the char-

acter of the problem. For example, if we are interested in the action of time-

dependent small perturbations, resulting from ¬‚uctuating external ¬elds,

on a system that is itself stationary, the basis functions will be station-

ary as well. On the other hand, if the basis functions are solutions of the

time-independent Schr¨dinger equation that contains moving parameters (as

o

nuclear positions), the basis functions are themselves time-dependent. The

equations for the mixing coe¬cients are di¬erent in the two cases.

In the following subsections we shall ¬rst consider direct evolution of

the wave function in a time-dependent ¬eld by integration on a grid (Sec-

tion 5.2.1), followed by a consideration of the evolution of the vector repre-

senting the wave function in Hilbert space. In the latter case we distinguish

two cases: the basis set itself is either time-independent (Section 5.2.2) or

time-dependent (Section 5.2.3).

5.2.1 Integration on a spatial grid

One way to obtain a solution is to describe the wave function on a spatial

grid, and integrate values on the grid in small time steps. This actually

works well for a single particle in up to three dimensions, where a grid of up

to 1283 can be used, and is even applicable for higher dimensionalities up

to six, but is not suitable if the quantum system contains several particles.

For the actual solution the wave function must be sampled on a grid.

Consider for simplicity the one-dimensional case, with grid points x0 =

0, x1 = δ, . . . , xn = nδ, . . . xL = Lδ. We assume that boundary conditions

for the wave function are given (either periodic, absorbing or re¬‚ecting) and

that the initial wave function Ψ(x, t) is given as well. We wish to solve

2 ‚2

‚Ψ(x, t) i

=’ ’ + V (x, t) Ψ. (5.6)

2m ‚x2

‚t

The simplest discretization of the second spatial derivative1 is

Ψn’1 ’ 2Ψn + Ψn+1

‚2Ψ

= , (5.7)

‚x2 δ2

yielding

2 2

‚Ψ(x, t) i

=’ + Vn Ψn ’ (Ψn’1 + Ψn+1 )

mδ 2 2mδ 2

‚t

1 The error in this three-point discretization is of order δ 2 . A ¬ve-point discretization with error

of order δ 4 is (1/12δ 2 )(’Ψn’2 + 16Ψn’1 ’ 30Ψn + 16Ψn+1 ’ Ψn+2 ). See Abramowitz and

Stegun (1965) for more discretizations of partial derivatives.

116 Dynamics of mixed quantum/classical systems

= ’i(HΨ)n , (5.8)

where Vn = V (xn , t) and

⎛ ⎞

b

a0

⎜b ⎟

a1 b

⎜ ⎟

⎜ ⎟

.. .. ..

⎜ ⎟

. . .

⎜ ⎟

H=⎜ ⎟. (5.9)

b an b

⎜ ⎟

⎜ ⎟

.. .. ..

⎜ ⎟

. . .

⎜ ⎟

⎝ b⎠

b aL’1

b aL

Note that we have absorbed into H. The matrix elements are

Vn

b=’

an = + , . (5.10)

mδ 2 2mδ 2

We seek the propagation after one time step „ :

Ψ(t + „ ) = e’iH„ Ψ(t), (5.11)

using an approximation that preserves the unitary character of the prop-

agator in order to obtain long-term stability of the algorithm. A popular

unitary approximation is the Crank“Nicholson scheme:2

UCN = (1 + iH„ /2)’1 (1 ’ iH„ /2), (5.12)

which is equivalent to

(1 + iH„ /2)Ψ(t + „ ) = (1 ’ iH„ /2)Ψ(t). (5.13)

Both sides of this equation approximate Ψ(t+„ /2), the left side by stepping

„ /2 backward from Ψ(t+„ ) and the right side by stepping „ /2 forward from

Ψ(t + „ ) (see Fig. 5.3). The evaluation of the ¬rst line requires knowledge

of Ψ(t + „ ), which makes the Crank“Nicholson step into an implicit scheme:

[2 + i„ an (t + „ )]Ψn (t + „ ) ’ ib„ [Ψn’1 (t + „ ) + Ψn+1 (t + „ )]

= [2 ’ i„ an (t)]Ψn (t) ’ ib„ [Ψn’1 (t) + Ψn+1 (t)]. (5.14)

Finding Ψ(t + „ ) requires the solution of a system of equations involving a