Lecture 3 - Quantum mechanics of few-state systems

by Leonid Pryadko

Schrodinger equation in continuum

ψ=.

eqs={I D[ψ[x,t],t]==-(1/2)D[ψ[x,t],{x,2}]+U[x] ψ[x,t],ψ[x,0]==(Subscript[ψ, 0])[x]}

Spread of a wavepacket

sol=NDSolve[Join[{ψ[-10,t]==0,ψ[20,t]==0},eqs/.{U[x]->0,(Subscript[ψ, 0])[x]->Exp[-(x-2)^2/2]}],ψ,{x,-10,20},{t,0,5}]

Manipulate[Plot[Evaluate[{Re[ψ[x,t]],Abs[ψ[x,t]]^2}/.sol[[1]]],{x,-10,10},PlotRange->All,PlotLabel->"t="<>ToString[t]],{t,0,10}]

Harmonic oscillator: time-dependent solution

sol=NDSolve[Join[{ψ[-10,t]==0,ψ[10,t]==0},eqs/.{U[x]->x^2/2,(Subscript[ψ, 0])[x]->Exp[-(x-2)^2/2]}],ψ,{x,-10,10},{t,0,10}];

Manipulate[Plot[Evaluate[{Re[ψ[x,i]],Abs[ψ[x,i]]^2}/.sol],{x,-10,10},PlotRange->All,PlotLabel->"t="<>ToString[i],PlotStyle->{Hue[0.0],Hue[0.4]}],{i,0,10}]

Double-well potential: time-dependent solution

F[x_]:=If[x<0,(x+2)^2/2,(x-2)^2/2]

Plot[F[x],{x,-4,4}]

sol=NDSolve[Join[{ψ[-10,t]==0,ψ[10,t]==0},eqs/.{U[x]->F[x],(Subscript[ψ, 0])[x]->Exp[-(x-2)^2/2]}],ψ,{x,-10,10},{t,0,60}];

Manipulate[Plot[Evaluate[{Re[ψ[x,t]],Abs[ψ[x,t]]^2}/.sol[[1]]],{x,-10,10},PlotRange->All,PlotLabel->"t="<>ToString[t],PlotRange->{-1,1},PlotStyle->{Hue[0.0],Hue[0.4]}],{t,0,60}]

Eigenvalue problems complicated

eqs1={EE ψ[x]==-(1/2)D[ψ[x],{x,2}]+U[x] ψ[x]}

sol1=NDSolve[Join[{ψ[-10]==0,ψ'[-10]==0.01},eqs1/.{U[x]->F[x],EE->%[[1]]}],ψ,{x,-10,10}];Plot[Re[ψ[x]]/.sol1,{x,-10,10}]

Need to dyscretize

Finite Difference solution

Let NN be the number of points per unit length

NN=25

This is the list of values

Table[A[i],{i,-10 NN-1,10 NN+1}];

This is kinetic energy in finite differences

KE=-(NN^2/2)ListConvolve[{1,-2,1},Table[A[i],{i,-10 NN-1,10 NN+1}]]/.{A[-10 NN-1]->0,A[10 NN+1]->0};

PE=Table[F[i/NN] A[i],{i,-10 NN,10 NN}];

vars=Table[ A[i],{i,-10 NN,10 NN}];

KE+PE

CoefficientArrays[KE+PE,vars]//Normal

matrix=CoefficientArrays[KE+PE,vars][[2]]//Normal;

Eigenvalues[matrix//N]//Sort

{vals,vecs}=Eigensystem[matrix//N]; vals

ord=Ordering[vals]

vals=vals[[ord]]; vecs=vecs[[ord]];

vals

Conjugate[vecs[[1]]].(matrix).vecs[[1]]

Conjugate[vecs[[2]]].(matrix).vecs[[2]]

{vals[[1]],vals[[2]]}

Two very close eigenvalues

v1=vecs[[1]];v2=vecs[[2]];

ListPlot[{v1,v2},PlotRange->All,Joined->{True,False}]

Single-spin dynamics

Solve Schroedinger equation using DSolve

eqs={I D[y1[t],t]==H11 y1[t]+H12 y2[t],I D[y2[t],t]==H21 y1[t]+H22 y2[t]}

DSolve[eqs,{y1[t],y2[t]},t]//Simplify

sol=DSolve[eqs/.{H11->1/2,H22->-1/2,H12->0,H21->0},{y1[t],y2[t]},t]

Average components of the spin: Introduce Pauli matrices

(Subscript[σ, x])={{0,1},{1,0}};(Subscript[σ, y])={{0,-I},{I,0}};(Subscript[σ, z])={{1,0},{0,-1}}; ψ={y1[t],y2[t]};

{(Subscript[σ, x])//MatrixForm,(Subscript[σ, y])//MatrixForm,(Subscript[σ, z])//MatrixForm}//MatrixForm

(Subscript[σ, x]).(Subscript[σ, x])+(Subscript[σ, y]).(Subscript[σ, y])+(Subscript[σ, z]).(Subscript[σ, z])//MatrixForm

(Subscript[S, x])=Conjugate[ψ].((Subscript[σ, x])/2).ψ;
(Subscript[S, y])=Conjugate[ψ].((Subscript[σ, y])/2).ψ;
(Subscript[S, z])=Conjugate[ψ].((Subscript[σ, z])/2).ψ;

{(Subscript[S, x]),(Subscript[S, y]),(Subscript[S, z])}/.sol[[1]]

Simplify[%,Assumptions->t>0]

Simplify[((Subscript[S, x])-I (Subscript[S, y]))/.sol[[1]],Assumptions->t>0]

Simplify[((Subscript[S, z]))/.sol[[1]],Assumptions->t>0]

Equation generator

It is a nuisance to have to specify equations explicitly all the time.  It would be more convenient to have Mathematica do it.  It is not immediately obvious how to do it:

eqs= (I D[ψ,t]-((Subscript[σ, z])/2).ψ)

DSolve[eqs,ψ,t]

Here is a function to help with it

DoEq[x_]:=Block[{i,y={}},
For[i=1,i≤Length[x],i++,y=Append[y,(x[[i]]==0)]];
y]

DoEq[eqs]

DSolve[%,ψ,t]

DoEq[x_,psi_,psi0_]:=Block[{i,y={}},
For[i=1,i≤Length[x],i++,y=Append[y,(x[[i]]==0)]];
For[i=1,i≤Length[ψ],i++,y=Append[y,(psi[[i]]==psi0[[i]])]]; 
y]

ψ/.t->0

eq1=DoEq[eqs,ψ/.t->0,{1,0}]

DSolve[eq1,ψ,t]

eqs=DoEq[I D[ψ,t]-(((Subscript[σ, z])/2)+2/5((Subscript[σ, x])/2)).ψ,ψ/.t->0,{1,1}/Sqrt[2]]

sol=DSolve[eqs,ψ,t][[1]]

Plot[(Subscript[S, x])/.sol,{t,0,10}]

Plot[Evaluate[{(Subscript[S, x]),(Subscript[S, y])}/.sol,{t,0,10}]]

ParametricPlot3D[Evaluate[{(Subscript[S, x]),(Subscript[S, y]),(Subscript[S, z])}/.sol],{t,0,10}]

You can now easily solve any spin dynamics problem.  Here is the magnetic field with the components (Subscript[B, x])=cos t, (Subscript[B, z])=1:

eqs=DoEq[I D[ψ,t]-(Subscript[σ, z]).ψ+ Cos[t]((Subscript[σ, x])/2).ψ,ψ/.t->0,{1,2}/Sqrt[5]]

sol=First[NDSolve[eqs,ψ,{t,0,10}]]

ParametricPlot3D[Evaluate[{(Subscript[S, x]),(Subscript[S, y]),(Subscript[S, z])}/.sol],{t,0,10}]