the+evolution+operator

toc =Introduction= We have seen in the lecture how to treat equilibrium quantum problems, using the density matrix of a system in contact with a thermal bath of inverse temperature //β//. In this Class Session, we examine the different problem of finding the time-evolution of a wave function //Ψ//(//t//) which starts from an initial state //Ψ// 0 at time //t//=0, and evolves through the Schrödinger equation. math i \frac{\partial}{\partial t} \Psi(t) = H \Psi(t) \quad\quad\text{(with $\hbar=1$)} math We study in this Class Session how to tackle this problem using the evolution operator.
 * Class Session 05 : The Evolution Operator **

A naive method: time discretization
Solving the Schrödinger evolution by remplacing the time derivative by a first order approximation math \frac{\partial}{\partial t} \Psi(t) \simeq \frac{\Psi(t+\Delta t) - \Psi(t)}{\Delta t} math is simple to implement but raises problems (it is difficult to ensure the conservation of probability).

A better method: the time evolution operator
It is better to start from the operator solution of the Schrödinger equation math \Psi(t) = \exp( - i t H) \Psi(0) math We observe that, putting formally //β// = //i t// (that is, working in //imaginary time//), one recovers the same evolution as that of the density matrix for an equilibrium system at inverse temperature //β//.

=Implementing the time evolution= We consider a particle of mass //m//=1, momentum //p// and position //x//, in a potential //V//(//x//). Its wave function //Ψ//(//x//,//t//) evolves with the Schrödinger equation of Hamiltonian math H=\frac{p^2}{2} + V(x) math where //p// 2 has the following representation in direct space math p^2 = -\frac{\partial^2}{\partial x^2}\:. math

The Trotter formula
As in the equilibrium approach, the idea is to decompose the time interval //t// in //N// steps of duration //Δt// = //t// / //N//. The full operator of evolution can be written (//without approximation//) as a product of //N// operators math \exp(- i t H) =\prod_{k=1}^N \exp( - i ~\Delta t~ H) math where each infinitesimal operator of evolution writes math \exp( - i ~\Delta t~ H) =\exp\Big[ - i ~\Delta t~ ( \frac{p^2}{2} + V(x)) \Big]\:. math In this expression, the potential part //V//(//x//) takes a simple form in real space, while the derivative part //p// 2 takes a simple form in the momentum space (that is, after a Fourier transform of the space coordinate). To separate both contributions, we use the Trotter formula math \exp\Big[ - i ~\Delta t~ ( \frac{p^2}{2} + V(x)) \Big] \simeq \exp\Big[- i V(x)~\Delta t/2 \Big] ~ \exp\Big[-i ~\Delta t\frac{ p^2}{2}\Big] ~\exp\Big[ - i V(x)~\Delta t~/2\Big]\:. math To which order in //Δt// is that formula valid? Check.

The algorithm
A possible scheme for the algorithm is thus:
 * Start from an initial wave function.
 * Discretise the time //t// in //N// steps of duration //Δt// = //t// / //N .//
 * At each time step, apply the Trotter formula by:
 * 1) in real space, multiplying by //e// −i //Δt V//(//x//) / 2 ,
 * 2) going to momentum space and multiply by //e// − //Δt// //p//^2 / 2 ,
 * 3) going back to real space, multiply by //e// −i //Δt V//(//x//) / 2.

Notes for the implementation
math \Psi(p) = \int_{-\infty}^{\infty} dx \exp(- i p x) \Psi(x) math math \Psi(x) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} dp \exp( i p x) \Psi(p) math In the algorithm, the space coordinate is also discretized and one uses the discrete Fourier transform.
 * The continuous space Fourier transforms write
 * An array may be used to encode the wave function: this simplifies the operations.
 * The discrete Fourier transform assumes the system to be periodic in space. This means, that the algorithms in the examples below actually simulate a wave function in a periodic potential landscape.

=Examples=

Oscillating in the harmonic potential
code format="python" from pylab import * steps = 400 space = linspace(-20,20,steps) xspace = pspace = space      # real space and momentum space lattices delta = space[1]-space[0] T=40                         #final time t_step=0.5 def fourierxp(phix): return array([integral(phix * exp(-1j * p * xspace)) for p in pspace]) def fourierpx(phip): return array([1./(2*pi)*integral(phip * exp(1j * x * pspace)) for x in xspace]) def integral(integrand): return sum(integrand)*delta def time_step(psi0,pot): psi0 = exp(-1j * pot * t_step / 2) * psi0 psi0 = fourierxp(psi0) psi0 = exp(-1j * pspace**2 * t_step / 2) * psi0 psi0 = fourierpx(psi0) return exp(-1j * pot * t_step / 2) * psi0 def zeit_ent(psi0, pot): evaluation = psi0 t = 0 while t < T:       print t        t += t_step evaluation = concatenate((evaluation,time_step(evaluation[-1],pot))) return evaluation pot=array([ [x**2/10 for x in xspace] ]) psi=array([ [2. * 1./sqrt(2*pi) * exp(-((x-1)/2)**2 / 2) for x in xspace] ]) psi = zeit_ent(psi, pot) ion line, = plot(xspace, abs(psi[0])**2) plot(xspace, pot[0]/4) axis([space[0],space[-1],-.1,1]) for n in range(len(psi)): line.set_ydata(abs(psi[n])**2) draw show code



We consider an harmonic (i.e. quadratic) potential with an inital Gaussian wave function.
 * What is your prediction for the time evolution of the probability density |//Ψ//(//x//,//t//)| 2 ?
 * Check by coding the time evolution of the wave function in python.

The tunnel effect



 * What is your prediction for the time evolution of the probability density |//Ψ//(//x//,//t//)| 2 ?
 * Check by coding the time evolution of the wave function in python.

Interferences in a Gaussian trap



 * What is your prediction for the time evolution of the probability density |//Ψ//(//x//,//t//)| 2 ?
 * Check by coding the time evolution of the wave function in python.

=Example animations=

Interferences in a Gaussian trap


url}?f=print| [Print this page]